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ABSTRACT 


This  report  contains  the  theoretical  background  material  for  a 
linearized  internal  gravity  wave  normal  mode  computational  scheme,  y 
implemented  at  the  Defence  Research  Establishment  Pacif ie^? Currents 
predicted  by  this  model  have  been  used  extensively  in  the  production  of 
simulated  synthetic  aperture  radar  images  of  internal  wave  signatures 
produced  by  moving  vessels.  A  full  derivation  is  given  for  the  steady  state 
internal  wave  fields  produced  by  a  distribution  of  fluid  volume  sources 
undergoing  common  uniform  horizontal  motion  in  a  finite  depth  ocean.] 
Features  include  a  full  normal  mode  expansion  of  the  steady  stajte', 
incorporating  both  the  radiating  internal  wave  field  and  the  localized 
disturbance  near  the  source.  Effects  of  background  shear  are  not  considered. 

-A  modif ication  of  the  standard  eigenfunction  expansion  leads  to  accelerated 
convergence  of  the  modal  sums.  Computation  of  all  relevant  fluid  dynamical 
fields  (2-component  velocity,  density,  pressure  and  vertical  particle 
displacement)  is  considered.  Both  sub-  and  super-critical  source  speeds  can 
be  handled.  An  extensive  discussion  of  the  effects  of  interior  regions  of 
high  evanescence  for  certain  wavenumber  regimes  shows  how  the  eigenfunction 
expansion  can  be  readily  modified  to  handle  such  regions  in  a 
straightforward  manner.  The  method  properly  accounts  for  trapping  of  energy 
generated  within  multiple  thermoclines  and  may  be  of  interest  in  certain 
underwater  acoustics  problems,  where  the  mathematics  are  similar.  Efficient 
computation  of  extensive  sets  of  field  data  is  discussed.  Examples  are 
presented  for  a  submerged  source  in  the  N.E.  Pacific  Ocean.  Also  included  is 
a  brief  comparison  with  measured  surface-ship-generated  internal  wave 
current  daca. 
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1,  INTRODUCTION 


Much  of  the  research  effort  of  the  Fluid  Dynamics  Group  at 
D.R.E.P.  has  been  directed  towards  the  remote  sensing  of  internal  wave  wakes 
using  various  devices,  including  Synthetic  Aperture  Radar  and  an  Infrared 
Scatterometer1 .  A  necessary  adjunct  for  prediction  of  the  performance  of 
such  devices  has  been  a  theoretical  model  to  predict  realistic  internal  wave 
current  fields.  These  current  fields  may  then  be  used  to  compute  simulated 
remote  sensing  imagery.  Applications  are  beyond  the  scope  of  this  report, 
which  describes  the  theoretical  background  to  the  internal  wave 
computational  model. 

The  density  of  the  Earth's  oceans  is  generally  non-uniform,  and 
especially  so  in  depth,  due  to  the  effects  of  gravity,  temperature,  and 
salinity.  This  stratification  of  density  allows  the  possibility  of  internal 
wave  propagation  due  to  buoyancy  effects.  In  a  stably  stratified  fluid,  the 
density  generally  increases  with  depth.  Thus  a  fluid  particle  which  is 
displaced  slightly  up  (down)  will  be  heavier  (lighter)  than  the  surrounding 
fluid,  and  so  be  subject  to  a  restoring  buoyant  force  towards  its 
equilibrium  position.  Hence  any  perturbation  of  a  stably  stratified  ambient 
fluid  will  generally  produce  internal  waves. 

Generating  mechanisms,  both  natural  and  artificial,  may  be  of 
many  types.  The  mechanisms  of  interest  here  are  effects  due  to  moving 
submerged  bodies,  which  can  lead  to  production  of  internal  waves  in  several 
ways.  One  such  mechanism  is  "wake  collapse"2’3,  which  involves  mixing  of  the 
density  by  the  source,  followed  by  collapse  of  the  density  in  the  wake  to 
the  ambient  stratified  values  after  passage  of  the  source.  A  second 
mechanism  is  the  so-called  "hull  effect"  which  involves  the  displacement  of 
fluid  particles  from  their  equilibrium  levels  due  the  action  of  the  hull. 
This  will  be  used  for  illustrative  purposes  later.  However,  the  actual 
details  of  the  source  model  have  little  direct  bearing  on  the  normal  mode 
expansion  which  forms  the  core  of  this  report. 
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Full  computation  of  an  internal  wave  field  produced  by  a 
realistic  submarine  would  be  an  extremely  difficult  task.  It  would  require 
solution  of  the  full  nonlinear  fluid  dynamical  equations,  incorporating  such 
a  variety  of  effects  as  the  detailed  shape  of  the  source,  variations  in 
ocean  depth,  ambient  currents,  ambient  surface  wind  waves,  effects  due  to 
the  Earth's  rotation  and  lateral  boundaries,  non-uniform  source  motions,  and 
so  on.  In  order  to  compute  model  internal  wave  fields  at  a  reasonable  cost, 
approximations  are  necessary. 

It  will  be  assumed  throughout  that  the  ocean  is  an 
incompressible,  inviscid,  vertically  stratified4  fluid.  Effects  due  to  the 
Earth's  curvature  and  rotation  will  be  neglected,  as  will  non-uniform  ocean 
depths  and  background  currents.  Thus  the  basic  ambient  state  is  a  quiescent 
finite  depth  ocean  of  infinite  horizontal  extent. 

The  source  model  for  the  hull  effect  is  taken  to  be  a 
distribution  of  fluid  volume  sources.  It  is  well-known5  that  in  potential 
flow  theory,  elementary  distributions  of  such  sources  in  a  uniform  stream 
can  produce  closed  streamsurfaces  that  resemble,  and  so  may  serve  as  models 
for,  realistic  hull  shapes.  In  a  linearized  model,  where  internal  wave 

a 

effects  are  assumed  small,  it  is  expected  that  the  basic  model  reasonably 
approximates  the  desired  rigid  hull,  there  being  minor  perturbations  due  to 
the  stratification.  For  purposes  of  illustrating  the  methods  discussed  in 
this  report,  the  only  source  to  be  used  is  a  simple  Rankine  ovoid5,6. 
However,  the  internal  wave  model  has  been  used  extensively  with  other  source 
distributions. 

Further  approximations  include  linearization  of  the  fluid 
dynamical  fields  about  the  ambient  state,  and  the  Boussinesq 
approximation7,8.  Such  approximations  are  used  extensively4,6,9,10  for 
studies  of  this  type  and  provide  an  efficient  computation  scheme  via  a 
normal  mode  expansion  technique.  Finally,  it  is  assumed  that  the  source  is 
in  uniform  horizontal  motion,  so  that  transient  effects  due  to  startup  or 
manoeuvering  of  the  source  are  not  implemented  in  the  model.  Such  transient 
effects  are  briefly  discussed  in  the  text,  however. 
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This  report  is  then  basically  concerned  with  numerical 
implementation  of  the  above-mentioned  normal  mode  expansion.  As  noted 
earlier,  actual  source  details  are  of  little  consequence  for  the  discussion, 
and  an  elementary  source  as  mentioned  above  is  used  for  purposes  of 
i 1  lustration. 

Section  2  presents  the  details  of  the  model  approximations,  and 
the  pertinent  linearized  fluid  dynamical  equations.  The  physical 
perturbation  quantities  of  interest  are  the  pressure,  density,  velocity  and 
vertical  particle  displacement  from  equilibrium.  It  is  shown  how  these  may 
all  be  expressed  in  terms  of  the  displacement  field,  which  in  effect  serves 
as  a  generalized  potential.  This  solution  of  these  basic  equations  in 
Fourier  space  is  also  discussed. 

Section  3  treats  the  solution  of  the  basic  differential  equations 
of  the  displacement  Fourier  transform.  The  approach  is  based  on  a  Green's 
function  technique.  Partial  fraction  expansion  of  the  Green's  function  using 
complex  variable  theory  allows  extraction  of  the  singular  part  in  closed 
form.  The  remaining  sum  then  has  accelerated  convergence  over  the  customary6 
direct  eigenfunction  expansion.  The  section  concludes  with  expressions  for 
the  general  Fourier  space  solutions  with  an  arbitrary  source.  These  could 
serve  as  a  departure  point  for  full  time-dependent  computations  for 
arbitrary  sources  in  general  motion. 

Section  4  is  concerned  with  specialization  to  a  source  which  is 
abruptly  "turned  on1'  at  t=0  and  is  in  steady  uniform  horizontal  translatory 
motion  thereafter.  This  leads  to  a  simplification  in  the  general  solutions. 
It  is  noted  that  the  resulting  expressions  are  non-singular  and  so  could  be 
used  for  time-dependent  computations  including  transient  effects  due  to 
startup. 


In  Section  5,  expressions  are  derived  for  the  ultimate  steady 
state  fields.  This  is  accomplished  by  referencing  coordinates  the  moving 
source,  and  considering  the  asymptotic  forms  of  the  fields  as  time 
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progresses.  The  physical  field  solutions  are  reduced  to  double  Fourier 
inverse  contour  integrals. 

Section  6  is  devoted  to  inversion  of  the  along-track  Fourier 
inversion  via  residue  theory.  Some  preliminary  discussions  concerning 
eigenvalues  are  included,  to  be  expanded  in  Section  7.  Notes  concerning 
general  source  terms  include  consideration  of  the  consequences  of  mass 
conservation.  The  section  concludes  with  the  main  results  of  the  report, 
which  is  a  complete  eigenfunction  expansion  of  the  fields,  expressed  as  a 
single  cross-track  inverse  Fourier  transform.  The  expansion  includes  not 
only  terms  representing  the  usual  far-field  radiated  internal  wave  wake  due 
to  real  eigenvalues,  but  also  contributions  from  imaginary  eigenvalues  which 
contain  near-field  information  toward  the  source  extremities. 

Numerical  details  and  techniques  are  discussed  at  length  in 
Section  7.  After  a  general  discussion  of  solution  methods  for  the 
differential  equation  (which  is  necessary  for  computation  of  eigenvalues), 
the  model  is  restricted  to  a  multilayer  type,  and  the  advantages  discussed. 
In  particular,  the  method  may  be  thought  of  as  suitably  joined  local  WKB 
approximations  to  the  eigenfunctions,  rather  than  a  stepped  approximation  to 
a  physical  profile.  Problems  commonly  arise  in  modelling  realistic  profiles 
which  may  contain  interior  layers  or  groups  of  layers  which  are  highly 
evanescent  in  certain  wavenumber  regimes.  The  problems  and  resolution  are 
illustrated  using  a  simple  symmetric  three-layer  model.  It  is  explained  how 
such  regions  can  cause  both  multiple  eigenvalues  and  exponential  overflow  on 
a  digital  computer.  The  exponential  overflow  is  easily  handled.  Simple 
methods  are  discussed  for  stable  computation  of  slightly  modified 
eigenfunctions  in  the  presence  of  such  multiple  zeros.  Basically  the 
standard  mathematical  eigenfunctions  are  modified  by  regrouping  terms  in  the 
eigenfunction  expansion.  This  is  accomplished  by  choosing  interior  reference 
depths  for  integration  of  the  eigenfunctions  in  a  manner  dependent  on  the 
profile  and  particular  associated  wavenumber.  The  resulting  modified 
eigenfunction  set  properly  accounts  for  channelling  of  energy  in  interior- 
ducts  and  can  accommodate,  for  example,  multiple  sources  on  opposite  sides 
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of  highly  evanescent  regions.  A  brief  comment  is  included  on  the  effect  of 
such  regions  on  "global  matrix"11  methods  used  in  acoustics,  where  the 
mathematics  is  similar  to  that  in  the  internal  gravity  wave  problem.  The 
section  concludes  with  a  discussion  for  setting  up  the  integrands  required 
for  Fourier  inversion.  Comments  for  an  efficient  computational  scheme  for 
computing  extensive  sets  of  field  points  are  included. 

Section  8  contains  examples  of  internal  wave  fields  computed 
using  the  methods  discussed.  A  slightly  smoothed  version  of  a  stability  (N  ) 
profile  computed  from  measured  summertime12  CTD  data  taken  at  Ocean  Station 
P  (50°N,145°W)  is  used  for  the  examples.  This  profile  contains  two  distinct 
thermoclines  which  gives  rise  to  some  of  the  problems  related  to  interior 
evanescence  mentioned  above.  Results  are  presented  for  two  source  speeds, 
one  of  which  results  in  purely  diverging  waves,  whereas  the  other  generates 
a  transverse  component.  Examples  are  also  given  with  the  source  travelling 
in  either  of  the  two  thermoclines. 

A  brief  comparison  with  measured  data  is  considered  in  Section  9. 
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2.  MODEL  AND  BASIC  EQUATIONS 


Relative  to  Cartesian  coordinates  (x,y,z),  the  model  under 
consideration  consists  of  an  inviscid,  adiabatic  fluid  of  infinite 
horizontal  (x,y)  extent  and  of  finite  depth,  d,  in  the  vertical  z-direction. 
Coordinate  z  measures  depth  positively  downv ‘d  from  a  mean  free  surface  at 
z=0,  in  the  direction  of  a  uniform  gravitational  field  g.  For  time  t<0,  the 
fluid  is  at  rest  and  the  density  pQ( z)  is  stratified  in  depth.  The 
stratification  is  assumed  stable,  so  ..tat  the  density  gradient  p'0{z )  is 
nowhere  negative.  Additionally,  it  is  nssumed  that  this  gradient  is  bounded 
within  the  fluid,  thereby  requiring  that  the  ambient  density  be  at  least 
continuous  for  0<z<d.  At  time  t=0,  a  distribution  S  of  volume  sources  is 
introduced  into  the  fluid,  which  results  in  the  generation  and  subsequent 
evolution  of  an  internal  gravity  wave  field.  This  report  is  concerned  with 
computation  of  the  associated  fluid-dynamical  fields  in  a  linear 
approximation. 

The  ambient  or  background  state  is  characterized  by  the  density 
p0(z)  and  pressure  p„(z)  which  are  related  by8  p;(z)  =  g  p0( z),  the  prime 
denoting  differentiation  with  respect  to  z.  For  t>0,  the  dynamical  fields  of 
interest  are  the  total  pressure  pT  a  p0(z)  +  p,  the  total  density  pT  a  p0  + 
p,  the  x  — ,  y- ,  and  z-components  of  velocity  denoted  respectively  by  u,  v 
and  w,  and  the  vertical  displacement  (positive  downwards)  {  of  fluid 
particles  from  their  equilibrium  positions.  In  a  linear  approximation,  the 
perturbation  fields  (p,p,u,v,w,f)  are  assumed  sufficiently  small  that  any 
terms  involving  products  of  two  or  more  of  these  quantities  may  be  neglected 
in  the  dynamical  equations.  The  resulting  linearized  equations  are8  : 


(1) 

Conservation  of  mass  (volume) 

:  V'U  =  S, 

(2.1a) 

(2) 

Invariance  of  density  : 

dp/d  t  +  p'0{  z)w  =  0, 

(2,1b) 

(3) 

Conservation  of  momentum  : 

p09u/9t  +  Vp  -  pgz  =  0,  and 

(2.1c) 

(4) 

Definition  of  displacement  : 

w  =  8{/flt. 

(2. Id) 
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In  these  equations,"  denotes  a  unit  vector  along  the  appropriate  coordinate. 
At  this  stage,  all  perturbation  fields  are  functions  of  the  four  variables 
(x,y,z,t) . 


The  velocity  and  density  perturbations  are  readily  eliminated 
from  Equations  (2.1),  leaving  a  pair  of  coupled  differential  equations 

3p/3z  =  -p0( z)  (  N2  +  32/8t2)  (2.2a) 

and 

p0(z)  d^/dt2dz  .  v2  p  +  p0(z)  as/at  (2.2b) 

for  the  pressure  and  displacement.  In  Equations  (2.2),  V2  3  82/3x2+82/3y2 
denotes  the  horizontal  Laplace  operator,  and  the  Brunt-Vaisala  (or 
intrinsic)  frequency  N  is  defined8  by 

N2(z)  »  g  p;(z)/p0(z).  (2.2c) 

This  fundamental  quantity  (N)  represents  the  local  natural  frequency  of 
oscillation  of  a  fluid  particle  slightly  displaced  from  Its  equilibrium 
position,  and  is  real  on  account  of  earlier  assumptions. 

For  present  purposes,  it  proves  convenient  to  eliminate  pressure 
from  (2.2),  leaving  a  single  equation 

*>„  (aY/et2)'  .  v2  (N2 .  a2/at2)  f  .  (Poas/et)-  (2.3) 

for  the  displacement.  The  coefficients  depend  only  on  z,  so  subsequent  work 
is  simplified  by  application  of  Fourier  transforms  In  x,  y,  and  t. 
Specifically,  the  triple  transform  of  a  function  f(x,z,t)  is  defined  by 


f (k,z,o)  b 


J. 

2  n 


//a2, 


Ik'S 


/  dt  e' 
o 


■i«t 


—  oo 


f(x,z,t), 


(2.4a) 
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and  the  inverse  operation  is  given  by 


00  ,  00—  1  £ 

f(x,z,t)  =  //  d2k  e'lk‘x  /  du  e1wt  ?(k,z,o).  (2.4b) 

4rr  -®  -  ®-ie 

where  x  =  xx+yy  and  k  s  kxx+kyy  denote  horizontal  position  and  wave-number 
vectors  respectively,  and  boldface  Roman  type  denotes  a  two-component 
vector.  The  wave  vector  k  is  initially  real,  whereas  the  frequency  w  has  a 
negative  imaginary  part  to  ensure  convergence  of  (2.4"'  'or  'reasonable' 
source  functions.  Then  ?  will  be  an  analytic  function  of  w  in  the  half-plane 
Im  u  ^  s,  and  Equation  (2.4b)  then  guarantees  quiescence  prior  to  t=0  in 
keeping  with  th8  imposed  initial  conditions. 

Q 

At  this  point  the  customary  Boussinesq  approximation  is  adopted. 
Thus  the  generally  small  terms  p'0(z) / p0(z)  are  neglected  in  the  dynamical 
equations,  except  when  multiplied  by  g  (cf  N2),  in  which  case  they  represent 
the  buoyancy  forces  which  allow  the  internal  waves  to  propagate.  Equation 
(2.3)  represents  the  general  case.  Under  the  Boussinesq  approximation  it 
simplifies  and  Fourier  transforms  to  the  simple  second-order  ordinary 
differential  equation 

?"  -  k2  [1  -  N2(z)/b2]  ?  .  S'/(io).  (2.6a) 

Complete  specification  of  the  problem  requires  boundary  conditions.  An 
appropriate  condition  at  the  flat  impermeable  bottom  z=d  is14  one  of  zero 
particle  displacement.  Under  the  Boussinesq  approximation,  attention  is 
focused  on  purely  internal  wave  motions,  which  will  have  generally 
negligible  amplitudes  at  the  free  surface.  The  appropriate  condition  is13 
therefore  one  of  vanishing  displacement  at  z=0  also.  The  boundary  conditions 

?  =  0,  (  Za0 ,  d  )  (2.5b) 

and  differential  equation  (2.5a)  form  a  boundary  value  problem  in  the 
variable  z  for  the  particle  displacement  transform  All  other  fields  may 
be  computed  from  the  solution  to  (2.5)  using  the  following  equations 


'SIP 


obtained  from  the  transformed  versions  of  Equations  (2.1)  and  (2.2)  : 


||PI 


Pressure 


Horizontal  velocity  : 
Vertical  velocity  : 
Buoyancy  force  : 


-i«  p0(z)  (  i»?'  -  S  )  /  k2, 
k  P  /  (»P0) » 
iw  and 

-Po  ^  ?• 


(2.6a) 

(2.6b) 

(2.6c) 

(2.6d) 


At  this  point,  note  that  it  follows  immediately  from  (2.6d)  and  (2.4b)  that 
in  the  physical  domain, 


g  p(x,z,t)  =  -p0(z)  N2(z)  t(x,z,t). 


(2.7) 


i.e.  the  buoyancy  is  directly  proportional  to  the  displacement  at  any  fixed 
depth,  and  will  not  be  given  further  consideration. 
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3.  GENERAL  SOLUTIONS  OF  THE  BASIC  EQUATIONS 

The  inhomogeneous  boundary  value  problem  (2.5)  can  be  solved  in  a 
standard  manner  using  a  Green's  function4,15  approach.  The  Green's  function 


appropriate  to  the  present  problem  satisfies  the  system 

[  5?/dz2  +  Q(k,c;z)  ]  G(k,c;z\z)  *  d(z-z'),  (3.1a) 

G(k,c;z\0)  =  G(k,c;z*  ,d)  =  0,  (3.1b) 

where 

Q(k,c;z)  =  N2(z)/c2  -  k2,  (3.2) 

and  the  horizontal  phase  speed  c  is  defined  by 

c2  a  a2  /  k2.  (3.3) 

Solutions  of  the  associated  homogeneous  boundary  value  problem 

P'(z)  +  Q(z)  0( z)  =  0  ;  0(0)  =  0(d)  =  0  (3.4) 


can  be  used  to  construct  the  Green's  function  as  will  be  shown  shortly.  For 
each  fixed  value  of  k2,  Equation  (3.4)  represents  a  standard16  Sturm- 
Liouville  boundary-value  problem  for  c-2  and  so  admits  a  denumerable 
infinity  of  eigenvalues  c“2(k)  and  associated  eigenfunctions  0n(k,z).  The 
eigenvalues  form  an  ascending  sequence  for  each  fixed  k2; 

c^2(k)  <  cj2(k)  <  ...;  c“2(k)-*»  (n->°°).  (3.5) 

The  eigenfunctions  are  orthogonal, 

It  N2(z)  0«(z)  0„(z)  dz  =  ^Ck)  (3.6a) 
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and  so  the  scaled  eigenfunctions 

$n(k,z)  =  <pn(k,z)  /  /?n(k)  (3.6b) 

constitute  an  orthonormal  set. 

The  Green's  function  defined  by  (3.1)  can  be  written  in  the  form 

G(k,c;z\z)  =  0o(k,c;z<)  ^d(k,c;z>)  /  A(k,c),  (3.7) 

where  <p0  and  #d  satisfy  surface  and  bottom  "initial-value"  problems 
specified  by  the  homogeneous  differential  equation  (3.4a)  together  with  the 
single-point  boundary  conditions 

0o(k,c;Q)  =  0,  0;(k,c;Q)  =  1,  (3.8a) 

$d(k,c;d)  »  0,  $d(k,c;d)  =  1,  (3.8b) 

and  where  z*.  =  min(z,z‘)  and  z>  -  max(z,z').  The  Wronskian 

A(k,c)  s  0o(k,c;z)  ^d(k,c;z)  -  0;(k,c;z)  0d(k,c;z)  (3.9) 

is  independent  of  z,  and  so  has  the  alternative  expressions 

A(k,c)  =  ^0(k»c*»d)  =  -0d(k’c;°)  (3.10) 

on  account  of  the  boundary  conditions  (3.8). 

2 

For  each  k  ,  G  is  a  meromorphic  function  of  c  with  simple  poles 
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at  zeros  of  the  Wronskian  A.  Consequently,  it  admits  a  partial  fraction 
expansion15 


G(k,c;z‘  ,z)  =  Gw(k;z'  ,z) 


+ 


-  ei(k) 

S  -5—“ — ? 

n=l  c  (k)  -  c 
n 


In  this  expression 


^0(k»c;z<)  ^d(k,c;z>) 
3  A(k,c)  /  a  c"2 


c=cn (k) * 


(3.11) 


G^kjz'.z)  »  sinh(kz<)  s inh[k(z>-d) ]  /  [  k  sinh(kd)  ]  (3.12) 

is  the  solution  of  (3.1)  when  c"2  =  0  (  or  in  the  unstratified  case  N  -  0  ). 
The  numerator  terms  within  the  brackets  of  (3.11)  are  proportional  to  the 
eigenfunctions  : 


0o(k,cn;z)  =  *B(k.z)  ^(k,cn;0)  /  ^(k,0),  (3.13a) 

9Jd(k,cn;z)  =  9>n(k>z)  ^d(k»cnl°)  1  0n(k>°)-  (3.13b) 

The  Wronskian  derivative  in  the  denominator  of  (3.11)  can  be  expressed  in 
terms  of  an  integral  as  follows17.  Let  denote  differentiation  with 
respect  to  a  parameter  (not  z)  to  be  made  definite  later.  Then  from  the 
differential  equation  (3.4)  and  its  parametric  derivative,  the  following 


results  may  be 

derived  ; 

(*d 

fo  -  <t> 0  0d)' 

\ 

3  Q0o0d* 

(3.14a) 

(*• 

kv 

=  Q0d^d* 

(3.14b) 

Integration  with  respect  to  z  of  (3.14b)  over  0<z<z‘  and  of  (3.14a)  over 
z'<z<d,  and  subsequent  addition  of  the  results  shows  that 


It  Q0o0ddz  a  A  ♦  *d(d)*;(d)-*#(d)to(d)  $o(0)fl,(0Wd(0)$;(0). 
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The  final  four  terms  all  vanish  due  to  the  boundary  conditions  (3.8), 
resulting  in  the  simple  integral  relationship 

A(k , c)  =  /J  Q(k,c;z)  30(k,c;z)  3d(k,c;z)  dz.  (3.15) 

If  the  unspecified  parameter  is  taken  to  be  c  2,  then  it  is  found  that 
SA/Sc-2  =  f*  N2  0O  <f>d  dz.  Substitution  into  (3.11)  and  subsequent 
application  of  (3.13)  and  (3.6)  leads  to  the  (modified)  eigenfunction 
expansion  : 


G(k,c;z',z)  =  Go^kjz'.z)  +  GR(k,c;z',z)  '  (3.16a) 

where 

GR(k,e;z\z)  =  S  cj(k)  $  (k,z)  $  (k,z‘)  /  [  c2(k)  -  c2] 

R  n-1  n  n  n  n  (3.16b) 

and  where  GM  is  given  in  closed  form  by  (3.12).  Application  of  a  direct 

<  A 

eigenfunction  expansion  '  solution  of  (3.1)  leads  to  an  expression 
identical  to  (3.16),  but  with  the  leading  term  GM  absent,  and  the  numerator 
factor  c *  replaced  by  c2c2,  i.e. 


00 

G  (k,c;z‘,z)  =  E  c2  c2(k)  ?n(k,z’)  $n(k,z)  /  [  c2(k)  -  c2] 
n»l 

Combination  of  this  result  with  (3.16)  leads  to  the  identity 

G  (k;z',z)  -  E  c2(k)  £  (k,z)  3 n(k,z')  . 

n-1  n  n  n  (3.17) 

The  left-hand  side  is  independent  of  the  profile  N  by  (3.12).  This  identity 
could  be  ■  sed  to  check  the  eigenfunctions  c  iputed  for  a  general  N  profile. 

The  advantages  of  (3.16)  over  the  standard  expansion  are  its  more 
rapid  convergence,  and  extraction  of  the  singular  behaviour  in  z.  The  term 
G„  alone  accounts  for  the  z - z ' )  term  in  (3.11).  Standard  asymptotic 
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results  for  fixed  k  and  n-*»  show  that 

C;2  ~  n27T2/f2(d)  and  (3.18a) 

0n ( k t, z )  ~  sin[  nfff(z)/f(d)]/VN(z),  (3.18b) 

where 

f(z)  *  /;  N(z ' )  dz1.  (3.18c) 

Thus  cn  aO(n-1),  $n=0(l),  <j>‘n  *  0  ( n ) ,  etc.,  and  the  convergence  rate  of 
(3.16)  is  proportional  to  n“*.  Consequently  up  to  two  z  or  z'  derivatives 
may  be  taken  'under  the  summation  sign*  in  (3.16)  while  retaining  absolute 
convergence  of  the  resulting  sum.  The  convergence  rate  of  the  standard 
expansion  is  n“2.  A  single  derivative  reduces  this  rate  to  conditional,  and 
a  second  derivative  reduces  the  convergence  to  purely  oscillatory,  as  is 
required  to  produce  the  5-singularity  of  (3.1).  This  singular  nature  is 
carried  in  the  isolated  term  of  (3.16),  a  fact  reflected  in  the  enhanced 
convergence  rate  of  the  associated  sum. 

In  terms  of  the  Green's  function,  the  desired  solution  to 
Equation  (2.5)  is 

?(k,z,u)  =  fl  G(k,c;z,z' )  S'(k,z',»)  /  (iw)  dz'.  (3.19) 

This  is  readily  derived  from  Equations  (3.1)  and  (2.5),  and  subsequent 
integration  by  parts  leads  to  the  expression  : 

?(k,z,o)  =  -/J  <3G(k;z,z' )/Sz‘  S(k.z,»)/  (i«)  dz'  (3.20) 

The  differentiated  sum  in  GR  of  (3.20)  has  an  n-3  convergence  rate  and  can 
be  differentiated  again,  as  required  for  evaluation  of  the  pressure  from 
(2.6a).  As  may  be  seen  from  Equations  (2.6),  the  remaining  fields  are  very 
simply  evaluated  from  the  pressure  and  displacement  fields.  In  the 
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evaluation  of  the  pressure,  there  arises  the  term 

fajdzdi'  =  kcosh(kz<)cosh[l((z>-d)3/sinh(kci)  -  <5(z-z‘); 

the  delta  function  removes  the  isolated  $  term  in  (2.6a). 

The  solutions  for  the  triple  Fourier  transforms  of  the  various 
fields  can  all  be  expressed  in  the  common  form  : 

^(k.z.M)  =■  /J  S(k,z',«)  Kf(k,o?;z,z')  dz‘,  (3.21) 

where  f  denotes  any  one  of  the  six  fields.  The  various  kernels  are 


Mk,o;z,z') 

•i/w, 

(3.22a) 

\ 

\  =  M.(k,c;z,z')  x 

Kw(k,«;z,z') 

1  l-l. 

(3.22b) 

Kp(k,«;z,z') 

2 

fiwpQ/k, 

(3.22c) 

Ku(k,w;z,z’) 

»  M2(k,c;z,z‘ )  x  ikx/k2 , 

and 

(3.22d) 

Kv(k,«;z,z‘) 

Mk  /k2. 

(3.22e) 

Auxiliary  function  is  just  9G/3z‘  and  M2  is  the  regular  part  of  the 
second  mixed  derivative  92G/9z9z';  i.e.  with  the  3-function  arising  from 
differentiation  of  GM  as  noted  above  excluded.  Function  M1  has  the  explicit 
expression 

CO 

M^k.cjz.z')  =  S  rn(k;z,z')  cj(k)  /  [c2(k)  -  c2]  +  rQ(k;z,z') 

n=1  (3.23a) 

where 

rn(k;z,z')  -  *n(k,z)  ^(k.z1),  (3.23b) 

and 

r„(k;z,z')  .  |sinh(kz>  cosh[k(z,-d)]l/sinh(k(i).(z<z'>  (3.23c) 
^cosh(kz’)sinh[k(z  -d)]'  ( z > z 1 ) 


-16- 


Similarly,  M2  is  given  by 


where 

and 


M2(k,c;z,z')  =  S  r;(k;z,z‘)  cj(k)  /  [c2(k)  -  c2]  +  r^(k;z,z* ) 

nsl  (3.24a) 

r^(k;z,z*)  ^  ^(k,z)  ^(k.z1),  (3.24b) 

r^(k;z,z‘)  s  k  cosh(kz)  cosh[k(z'-d)]  /  sinh(kd).  (3.24c) 


Equations  (3.2i)  through  (3.24)  represent  the  general  transform 
space  solutions  for  the  various  fields  of  interest.  These  equations  could  be 
used  for  fully  time-dependent  computations  involving  arbitrary  source 
distributions  undergoing  arbitrary  motions  under  the  linear  and  Boussinesq 
approximations,  and  are  given  for  future  reference  purposes.  Later,  they 
will  be  specialized  somewhat  for  the  purposes  of  this  report. 
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4.  UNIFORMLY  TRANSLATING  SOURCE:  FREQUENCY  INVERSION 

At  this  point,  the  model  is  specialized  to  forcing  by  a 
distribution  of  fluid  volume  sources  over  some  finite  planar  region  at  depth 
z9,  in  uniform  common  translation  at  speed  c0  parallel  to  the  x-axis,  moving 
in  the  direction  of  increasing  x  and  "switched  on"  at  time  t=0.  Thus  the 
source  function  has  the  form 

S(x,z,t)  -  <5(z-zt)  $0(x-c0t,y)  H(t),  (4.1) 

where  co>0.  It  should  be  noted  that  several  such  sources  can  be  readily 
accommodated  by  superposition  of  the  individual  solutions.  In  a  refined 
scheme,  more  general  vertical  source  distributions  could  be  incorporated  by 
retaining  the  integration  over  the  Green's  function,  but  this  will  not  be 
considered  here.  The  triple  transform  of  the  speeial  source  (4.1)  is 

5(k,z,w)  =»  -  i  SQ(k)  <5(z-z$)  /  (w-cQkx),  Im{u-cokx}  <  0,  (4.2a) 


where 

00 

S  (k)  »  (2H)'1  U  S  (x)  exp(ik-x)  d2»  (4.2b) 

-00 

is  a  double  spatial  Fourier  transform.  The  general  expression  (3.21) 
immediately  simplifies  to 

F(k,z,«)  =  -i  S0(k)  Kf(k,w;z,z8)  /  (a>-c0kx).  (4.3) 

Recovery  of  the  solutions  in  the  physical  domain  then  involves 
inversion  of  the  general  expression  (3.21)  according  to  (2.4b).  For  the 
uniformly  translating  source,  it  proves  convenient  to  transform  to 
coordinates  x,  a  x  +  c0t  x  moving  with  the  source.  The  physical-space 
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solutions  then  have  the  general  form 


f(xrz.t)  .  J7  S0(k)  Lf(k,t;2,2s)  e"1*-*  d2k 


(4.4a) 


where 


4  «A|/  i  y  w_ v»rtixy  y  v 

Lf(k,t;2,2  )  ■  jij  /  Kf(k,ffli2,2s)  k -  do).  (4.4b) 

-co-15  0  X 


exp  i  (co-c  k  )t 


The  first  step  is  evaluation  of  the  frequency  inversion  integral  (4.4b)  for 
the  five  kernels  (3.22).  With  the  special  choice  (4.2a)  of  source,  it  proves 
possible  to  evaluate  these  integrals  using  residue  theory.  The  source  term 
introduces  a  simple  pole  at  tf=c0kx.  The  displacement  (3.22a)  contains  an 
additional  pole  at  w=0.  Moreover,  all  field  transforms  have  poles  where  a2  = 
k2  c2,  as  is  evident  from  Equations  (3.23a)  and  (3.3).  At  this  point  k  is 
real,  and  so  k2  is  merely  a  positive  parameter  as  far  as  a  inversion  is 

1  A 

concerned.  Standard  Sturm-Liouville  theory  then  shows  that  the  eigenvalues 
c~2  are  real  and  therefore  that  the  frequency  poles  ±kcn  are  on  the  real  o- 
axis.  There  are  no  other  poles  in  the  o-plane.  All  eigenfunctions  are 
independent  of  a.  Thus  all  poles  are  simple  and  appear  explicitly  as 
denominator  factors  of  the  form  (w±cjp).  The  frequency  inversion  contour 
passes  below  the  real  axis  and  so  all  fields  vanish  for  t<0  as  may  be  seen 
by  closing  the  inversion  contour  in  the  lower  half  plane.  For  t>0,  the 
contour  may  be  closed  up,  thereby  picking  up  contributions  from  the  poles  on 
the  real  axis.  Evaluation  by  residue  theory  is  simple  if  tedious,  and  leads 
to  the  following  results. 

(a)  Particle  displacement: 

Lf(k,t;z,z$)  =  i  rQ(k;z,zs)  *(©Q,t)  - 

t  S  e2(k)  rn(k[2,2s)  {[*(«-  t)]/Z  -  *(»  ,t)J 

n=l  r*  \ 


-19- 


(b)  Vertical  velocity: 


Lw(k,t;z,z$)  =  -  r0(k;z,z$)  + 

S  k  cj(k)  rn(k;z,zs)  [»(«",t)-t(eJ,t)]/Z 


(4.5b) 


(c)  Pressure: 


Lp(k,t;z,zs)  =  iu0/>0(z)r’(k;z,zs)/k  - 

ip0(z)  S  c*(k)  r‘(k;z,z$)  [♦(»jj,t)+*(oJJ,t)]/2 
n=l 


(d)  Horizontal  velocity: 


Lu(k,t;z,zs)  =  i  k  I^(k;z,z$)  /  k  • 

1  i  E  cj(k)  r'(k;z,z  )  [t(e',t)-*(oJ,t)]  /  2 
n=l 


In  these  equations,  note  the  auxiliary  definitions 


where 

and 


wo  3  cokx«  ®n  3  wo  ±  ®n» 

wn  3  k  c«(k). 

*(a,t)  s  [l-exp(-iot)]/a. 


(4.5c) 


(4.5d) 

(4.6a,b) 

(4.7) 

(4.8) 


✓S 

Also,  (4.5d)  is  a  vector  form,  and  k  is  a  unit  vector.  Equations  (4.4)  and 
(4.5)  give  the  solutions  in  physical  space  for  all  time  for  the  particular 
source  (4.1),  and  can,  in  principal,  be  used  to  compute  the  time-evolution 
of  the  linearized  fluid-dynamical  fields,  i.e.  including  transient  effects. 
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5.  EXTRACTION  OF  STEADY-STATE  FIELDS 

Of  greater  interest  for  present  purposes  is  the  ultimate  steady- 
state  situation  viewed  in  the  (uniformly  translating)  source  reference 
frame.  Specifically,  the  steady-state  form  f°°  for  field  f  is  defined  by  way 
of 

f"(x,y,z)  s  lim  f(x1sy,z,t)  =  lim  f (x+c  t,y,z,t)  (5.1) 

t-K»  ■  t-w>  0 

It  should  be  noted  that  for  finite  t>0,  the  various  integrands  of  Equations 
(4.5)  are  well-defined,  for  i(cE,t)  is  an  entire  function  of  both  arguments. 
As  t -->«>,  the  increasingly  oscillatory  nature  of  the  exponential  part  in  each 
term  ¥[>p(k)  ,t] ,  at  all  k-points  ?iear  which  0p(k)*O,  effectively  removes 
these  contributions  from  the  integral  (e.g.  Riemann-Lebesgue  lemma),  thereby 
replacing  each  such  term  i  by  cj'^k)  near  such  points.  However,  k-points 
near  real  roots  k0  of  «p(k)=0  present  problems.  These  lead  to  poles  in  the 
resulting  steady-state  integrands  on  the  respective  axes  of  integration.  The 
attendant  ambiguities  can  be  circumvented  as  follows8.  Near  such  points  the 
ensemble  of  integration  may  be  shifted  into  the  respective  complex  planes,  ( 
i.e.  k0+5k,  5k  complex  )  for  any  finite  t  due  to  the  analyticity  of  the 
integrands.  On  the  displaced  contour  near  such  a  zero, 

*[»p(k),t]  a  [l-exp(-itA)]/A  (5.2a) 

where 

A  a  <5k-?k  «p(k0).  (5.2b) 

This  deformation  is  then  to  be  done  in  such  a  way  that  Im{A}<0  everywhere  on 
the  deformed  portions;  the  corresponding  exponential  term  of  }  decays  in 
time,  and  the  resulting  term  A  1  is  non-singular  on  the  deformed  contour  in 
the  steady-state  limit.  If  there  are  points  k0  on  the  real  axes  at  which 
also  Vk  o>p  =  0,  a  higher  approximation  must  be  used;  this  possibility  will 
be  considered  further  later. The  relevant  cup  are  given  in  (4.6),  and  have  k- 
gradients: 
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7k  a0  a(co,0)  ;  (5.3a) 

Vk  «;  3  (  c0  ±  kxcgn(k)/k,  ±ky  cgn(k)  /  k  ) 

■  (  [<5  -  cn(k0)cgn(k#)3/c0,  ±koycgn(k0)  /  k0  ).  (5.3b) 

The  final  expression  on  the  right-hand  side  is  evaluated  at  the  appropriate 
k0,  and  the  horizontal  group  speed  is  defined  as 

cgn(k)  =  d«i;(k)/dk.  (5.3c) 

It  will  be  shown  later  that  | cgn | < |cn|< |c0|  for  poles  on  the  real  kx-axis, 
so  that  (c^-cncgn)o^0.  The  possibility  of  equality  is  ignored  here  for  now, 
and  only  the  kx-contour  is  deformed,  leaving  ky  as  a  real  parameter  in  the 
complex  kx-integrals.  The  deformation  criterion  (comment  following  (5.2b)  ) 
is  then  satisfied  by  deforming  the  kx-contour  below  each  resulting  steady 
state  pole  on  the  real  kx-axis. 

The  resulting  expressions  for  the  steady-state  physical  fields 
can  be  written  in  the  general  form 


f*(x,z)  3  (Htt)"1  / 


4 


s0(k) 


Lf(k,z) 


e~ix*k  dk 


dk. 


(5.4) 


Here  the  ky-integral  is  confined  to  the  real  axis,  and  the  kx-integral  is 
taken  along  a  contour  K  initially  consisting  of  the  real  kx-axis,  indented 
below  any  real  poles.  With  dependence  of  the  various  wp  and  T  on  their 
respective  arguments  (k)  and  (k;z,z8)  suppressed  in  the  notation,  the 
kernels  have  the  expressions  : 

(a)  Particle  displacement: 

L>.z)  ■  V°o  -  1  *  rn  k2  '  t  ®0(  "o  ‘  “n  >  3- 
'  n=l 


(5.5a) 
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(b)  Vertical  velocity: 


Lw(k,z)  =  -  ro 


00 
k  E 
n=l 


r  k2  c*(k)  /  (  u2  -  a2  ). 
n  o'1  '  v  o  n  ' 


(5.5b) 


(c)  Pressure: 


>>•*)  ■  'Vori/k2  - 


ip^ca 
~0  0 


E 

n=l 


r  cj(k)  /  (  «2  -  <1  ). 

n  n'-  ’  v  o  n  ’ 


(5.5c) 


(d)  Horizontal  velocity: 


Lu(k,z)  .  i  k  r;/  k  - 


CO 

ik  E  T  k  c?(k)  /  (  <a2  -  w2  ). 
„  ,  n  nv  1  v  '0  n  ' 
n=  1 


(5.5d) 


Note  that  in  the  infinite  sums  for  w,  p  and  uh,  the  respective  convergence 
rates  in  n,  for  each  fixed  k  in  the  integrands,  are  0(n-3),  0(n~2),  Q(n-3). 
The  w0  term  in  the  infinite  sum  for  the  displacement  field,  which  may  be 
isolated  by  a  partial  fraction  expansion  of  the  denominator,  is  independent 
of  n,  and  by  (3.17)  this  isolated  term  could  be  evaluated  explicitly  to 
cancel  the  F0  term  outside  the  sum.  However  the  remaining  infinite  sum  would 
have  only  0(n-1)  convergence,  in  contrast  to  the  G(n-3)  behaviour  of  (5.5a). 


Finally,  note  from  (5.5c)  and  (5.5d)  that  Lp  a  c0  p0(z)  lu;  these 
expressions  are  essentially  partial  fraction  expansions  of  the  respective 
fields,  and  the  numerators  have  the  same  values  at  the  poles  <d0,  w*.  Thus 
the  steady-state  pressure  is  proportional  to  the  x-velocity  component  at  any 
depth.  This  may  also  be  seen  from  the  dynamical  equations  upon 
transformation  to  coordinates  translating  with  the  source  and  looking  for 
solutions  independent  of  time.  The  pressure  field  will  therefore  be  excluded 
from  further  consideration. 
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6.  EVALUATION  OF  k„-INTEGRALS 
(a)  Preliminaries 

As  in  the  case  of  the  frequency  integrals  earlier,  the  kx- 
integrals  (5.4)  arising  in  computation  of  the  steady  state  fields  are 
amenable  to  evaluation  by  residue  theory. 

There  are  poles  due  to  zeros  of  the  sinh(kd)  (k*0)  terms  in  the 
denominators  of  F0  and  There  are  also  poles  at  kx*0  in  the 
displacement,  and  a  pole  at  k=0  (  and  therefore  at  kx  »  ±i|kyl  Jin  the 
horizontal  velocity  fields.  These  poles  give  rise  to  near-field  disturbances 
and  are  independent  of  the  N2-profile  details.  Finally,  there  are  poles 
where  w2  =  <a2,  which  depend  on  the  details  of  the  profile  and  represent  the 
radiating  internal  wave  wake. 

The  coefficient  Q  of  £  in  the  differential  equation  (2.5a)  is 
[N2(z)/c2-k2] .  This  coefficient  must  be  positive  somewhere  on  (0,d)  in  order 
that  the  solution  not  be  of  a  purely  exponential  nature.  Thus  it  is  evident 
that  for  k2-*°°,  c2->0+.  It  has  already  been  seen  (in  the  discussion  for  g“) 
that  if  |c2|->«,  then  k2->  (nir/d)2.  Moreover,  when  k2<0,  the  d.e.  is  always 
oscillatory  for  c2>0,  and  also  for  certain  values  of  c2<0.  This  leads  to  the 
basic  rules  that  (k2>0)  =>  (c2>0),  and  conversely  (c2<0)  =»  (k2<0).  The  basic 
behaviour  of  the  eigenvalue  curves  c2  »  c2(k)  in  the  (c2,k2)  plane  is 
sketched  in  Figure  1  (lighter  curves  with  negative  slope  ).  These  curves 
were  computed  for  a  constant-N  ocean,  and  are  included  purely  for 
illustrative  purposes. 

2  2 

The  poles  &>n  =  u0  can  be  located  as  follows.  These  must  satisfy 
the  dual  constraints  c2  ■  c2k2/k2  (which  is  a  restatement  of  the  equation  u2 
a  u2  )  and  c2  =>  c2  (i.e.  they  must  be  eigenvalues).  For  a  fixed  ky,  the 
first  constraint  is  equivalent  to  a  hyperbola 

(c2-c2)k2  .  c2k2  (6.1) 
2  2 

in  the  (c  ,k  )  plane.  A  sample  of  such  a  hyperbola  is  also  sketched  as  the 
darker  curve  (positive  slope)  in  Figure  1.  The  intersections  of  the 


Figure  1.  Sample  dependence  of  five  eigencurves  cn  on  k  is  shown  by  the 
lighter  curves.  The  darker  curve  depicts  a  constraint  hyperbola  of  the  form 
(cj-c2)k2  =  c2k2. 

hyperbola  with  the  eigenvalue  curves  give  two  discrete  families  of  poles 
satisfying  the  constraint,  one  with  k2>0  and  the  other  with  k2<0.  Both 
families  have  a  horizontal  asymptote  c2=c2.  Poles  with  k2>0  have  0<e2<c2, 
whereas  those  with  k2<0  have  c2>c2.  In  both  cases  c2>0,  so  that  k2>0 
corresponds  to  k2>0,  i.e.  to  poles  on  the  real  kx-axis.  Similarly,  the  poles 
with  k2<0  correspond  to  poles  on  the  imaginary  kx~axis. 

In  order  to  evaluate  the  kx- 1 ntegr a  1  via  contour  integration  it 
is  necessary  to  compute  the  residue  at  each  pole  l/wn.  If  these  poles  are 
simple,  the  residue  calculation  involves  the  reciprocal  of  the  derivative 
dwn/dk  ,  evaluated  at  the  pole.  Since  the  modal  frequencies  depend  on  k  only 
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via  k2,  it  is  found  that 

da;/dkx  =  c0±(kx/6;n)dM2/dk2. 


(6.2a) 


Evaluation  at  «”=0  gives 

(d»*/dkx)0  »  [c2-(du2/dk2)0]/c0,  (6.2b) 

the  subscript  "o"  denoting  evaluation  at  the  pole.  Note  that  the  term 
dca2/dk2  can  also  be  written  as  cncgn  as  in  (5.3).  An  integral  expression  for 
this  term  can  be  obtained  as  follows.  From  (3.15)  and  (3.13)  it  follows 
immediately  that 

ft  Q(k,c;z)  $2(k,z)  dz  =  0,  (6.3) 

since  A(k,c)  *  0  along  any  eigencurve.  Note  that  the  derivative  in  (6.3)  is 
total,  since  there  is  just  one  free  parameter  describing  an  eigencurve.  Let 
the  parameter  be  k2,  then 

ft  [N2/a2  -  1  -  (k2N2/Qj)du2/dk2]  <f>2n  dz  -  0. 

This  is  readily  rearranged,  with  the  help  of  normalization  definition 
(3.6a),  to 


d«£/dk2  =  4  [1  -  1$  Jt  dz],  (6.4) 

The  left-hand  side  of  (6.2)  may  therefore  be  expressed  as 

(d«*/dkx)0  =  (c0/kj)[k2  .  c2  kj„  Jt  dz],  (6.5) 

The  term  in  square  brackets  is  strictly  positive  for  both  real  and  imaginary 
poles,  except  if  ky  =0,  and  kx=0  is  also  a  pole.  Evidently,  for  ky/0,  the 
poles  in  kx  are  all  simple.  Moreover,  even  if  ky=Q,  the  poles  are  still 
simple  so  long  as  the  corresponding  value  of  kxo  is  non-zero,  or  if  the 
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limiting  value  of  the  ratio  ky/kn  is  non-zero.  With  reference  to  Figure  1, 
note  that  as  ky-»0,  the  hyperbolic  curves  degenerate  to  the  lines  k2=0, 
c2=c2.  Clearly  a  double  pole  is  possible  ,.i  at  most  one  point,  and  then  only 
if  e2(0)  =  c2  for  some  eigencurve.  This  situation  could  be  resolved  by 
deforming  the  ky-contour  near  ky=0  also,  but  this  will  not  be  considered 
further  here.  Related  problems  are  removed  by  any  realistic  source  term  (see 
(6.11)  below).  Note  that  for  poles  with  k2>0,  the  right-hand  side  of  (6.4) 
is  then  positive,  as  mentioned  in  conjunction  with  (5.3c). 

From  the  definition  of  the  phase  speed,  Eq.(6.4)  can  be  expressed 
alternatively  as 

dej/dk2  *  -cj  fi  fn  dz.  (6.6) 

This  shows  that  each  curve  c2  is  a  strictly  decreasing  function  of  k2,  a 
standard7  result,  and  is  illustrated  by  the  sample  curves  of  Figure  1. 

For  the  real  roots  (k2>0),  the  constraint  c2k2  *  k2c2  may  be 
written  in  polar  form  by  defining  kx  *  k  cos 0,  ky  =  k  sin 9.  The  nth 
eigencurve  is  then  specified  in  the  form  (k,0n(k)).  The  constraint  reduces 
to 


cos2[0n(k)]  -  c2(k)/c2.  (6.7) 

The  right-hand  side  is  a  well-defined  function  of  k,  and  as  k->»,  it 
approaches  0,  and  so  6n-*±i\/Z;  the  curves  have  a  vertical  asymptote.  If  it 
happens  that  cn(0)  <  c0,  then  the  curve  passes  through  the  origin  in  the 
(kx,ky)-plane,  at  angle  f?n(0).  Such  curves  are  termed  super-Froude6  and  are 
infinite  in  number.  They  contain  only  diverging  waves.  Alternatively,  if 
cn(0)>co*  then  there  ^  sonie  finite  value  of  k  at  which  0n(O)=O.  The  slopes 
of  the  cirves  are  given  by 


dk/dfl  =  -cosin(0)/c^. 


(6.8) 
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Thi$  second  type  of  curve  is  termed  sub-Froude  and  each  such  curve  cuts  the 
kx-axis  with  vertical  slope  as  is  evident  from  (6.8)  -  they  give  rise  to 
modes  containing  transverse  waves.  There  are  at  most  a  finite  number  of 
modes  of  this  type.  If  the  source  speed  e0  is  sufficiently  high  that  c0  ^ 
0,(0),  then  all  modes  are  super-Froude. 

Similarly  the  constraint  for  the  imaginary  modes  may  be  rewritten 
as  follows.  Define  ks  »  ikcosh/i,  ky  =  ksinh/t  where  k^O.  Then 

cosh  2fi  =  cj(k)/cj.  (6.9) 

It  is  readily  seen  that  for  any  super-Froude  mode,  the  corresponding  /i  (and 
consequently  kx)  is  bounded  away  from  zero.  The  associated  wake  contribution 
thus  decays  exponentially  away  from  the  source  with  a  factor  exp(-|xkx(Q) | ) . 
Conversely,  any  sub-Froude  mode  has  /i-values  down  to  zero  and  the  associated 
contribution  from  the  kx-integrals  might  be  expected  to  decay  only 
algebraically  with  distance  from  the  source. 
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(b)  Source  notes 

The  source  (4.1)  is  represented  by  a  horizontally  bounded 
distribution  $0  of  sources  on  the  plane  z=z3.  Let  these  sources  lie  between 
xmln  and  xm<3X.  Then  in  general,  the  Fourier  transform  S0  will  possess 
exponential  growth  for  |Xm(kx)|-*».  Closure  of  the  various  kx-contours  will 
thus  be  possible  only  for  x  >  xm0X  and  x  <  xmln.  Thus  the  result  of  this 
contour  deformation  will  describe  the  far-field  propagating  internal  wave- 
field,  and  the  near-field  both  ahead  of  and  behind  the  extremes  of  the 
source  distribution.  Computations  for  xrain  <  x  <  xna!{  would  require 
alternative  computational  methods,  and  will  not  be  considered  further  here. 

Since  the  source  distribution  S0  is  real,  its  transform  (4.2b) 
satisfies  the  property 

C§(k)3*  -  s(-k*),  (6.10) 

the  asterisk  denoting  complex  conjugation. 

A  realistic  source  model  must  conserve  mass.  This  requires  that 
the  integral  of  V*u  taken  over  a  volume  completely  enclosing  the  source 
distribution  S  must  vanish.  For  the  source  (4.1),  this  requirement  reduces 
to  //r<»  S0(x)  d2x  =  0,  which  is  equivalent  to 

§(0)  »  0.  (6.11) 

This  constraint  will  be  assumed  in  what  follows.  It  alleviates  potential 
singular  behaviour  near  k»0. 

(c)  Evaluation  of  the  integrals 

It  has  been  seen  above  that  the  nth  term  of  each  infinite  sum  in 
Equations  (5.5)  contains  two  real  and  two  imaginary  poles  in  kx.  Denote 
these  as 


^nCky)’  ±il/n(ky)* 


(6.12) 
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where  nn,  vn  Z  0  by  definition.  As  noted  earlier,  there  are  additional  pole< 
at  zeros  of  sinh(kd)  (k^O).  These  lie  at  k2  =  -n27i2/d2  s  -y2,  and  so  in  the 
kx-plane  at 


k*  ■  -«5(ky)  ■  -  (  K?  ♦  ll  ) .  (6-13) 

where  by  definition,  dn  ^  0.  Thus  there  are  poles  at  kx  =  ±i<5n(ky).  The  plus 
(minus)  sign  thus  refers  to  poles  in  the  upper  (lower)  half  kx-plane. 
Finally  there  is  a  pole  at  kx*0  in  the  displacement  integral,  and  poles  at 
kx  =  ±i|kyl  in  the  horizontal  velocity  integrals. 

The  details  of  the  residue  calculation  then  form  a  standard  but 
tedious  exercise  in  algebra.  The  final  results,  which  represent  the  working 
equations  for  a  computational  scheme,  are  presented  below.  With  reference  to 
Equation  (5.4),  it  is  evident  that  all  fields  can  be  expressed  in  the  form 
of  a  single  inverse  Fourier  transform  in  ky. 

f“(x,z)  =  /“«  f(x,ky,z)  exp(-iyky)  dky,  (6.14a) 


where 


00 

f(x,ky,z)  S  (2 n)'1  Kf  SQ(kx,ky)  L”(kx,ky,z)  e~1xk*  dkx.  (6.14b) 


The  results  of  evaluation  of  the  kx-integral  (6.14b)  by  residue  theory 
(closing  up  if  x<xmin  and  down  if  x>xra<3X  as  necessitated  by  an  above  note) 
may  bo  summarized  as  : 

00 

f(x,k  ,Z)  a  E  ^fn(kv,Z,Zs) 
y  n=0  Tn  y  s 

n=l  Vyz,2s) 


af  V“i<5n’V  exP('x<5n)  " 
af  V^VV  exP('XI/n) 


(  x  >  x 


max 


) 


(6.15a) 
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f(x,ky,z)  =  Df(ky,z,z$)  SQ(0,ky)  + 

*  Hfn^y*z’zs)  af  V+i  W  exp^x<V  ’ 

n=0 


n=l 


f,  {  WV2'^  af  So(*lE'r>-ky)  e*P(«"n)  * 

iVyz,zs^  bf  SQ(+/tn,ky)  exp(-1x/tn)  + 

bf  V“VV  exp(+ixV  ) 

Efn(ky’z’zs)  V0,y  }* 


<  X  <  xmin  ) 


In  these  equations  the  a's  and  b's  are  constants  defined  by 

±  i  i  i 

au  ■  ±1,  av  =  +i,  aw  =  +1,  a^  =  ±1,  and 

b~  =  +1,  b~  S  ±1,  b*  a  ±\ ,  b^  »  +1. 

The  coefficients  H  (  Hyperbolic,  for  want  of  a  better  identifying 
arise  from  the  poles  in  the  sinh(kd)  term,  and  are  defined  by 


and 


lWy2,2s)  ■  1 

Vy2'2s>  ■  kyMn 


■  1 

W2'^  B  1/(W 


£ncos(7nz)cos(Vs)  /  d, 


7nsin(7nz)c°s(Vs)  7  (fy1) 


In  these  expressions, 


(6.15b) 

(6.16a) 
(6.16b) 
term  ) 

(6.17a) 

(6.17b) 

(6.17c) 
(6. 17d) 


en  a  (1/2,1)  (n=0,n>0). 


(6 . 17e) 
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The  coefficients  I  arising  from  the  poles  on  the  imaginary  kx-axis,  and  R 
arising  from  poles  on  the  real  axis,  due  to  the  terms  «n  of  (5.5),  share  the 
common  functional  form 


and 


Vyz’2,>  ‘ 1 

Gvn<yz-zs>  ■yiyi 


x 


cn  r;^-z’zs) 

2(c l  -  du^/dk2)  ' 


Gwn<ky-Z-2S>  *  lklcn 
G^.n(Ry,z,2s)  3  1 


Vn  rn(k;z’zs> 
2(c2  -  do2/dk2)  ' 


These  are  to  be  evaluated  with 


k„  =  *1vn{ky),  k  =  .Hul-W'2,  (6.1), 


n 


and 


k„  =  .  nM,  k  =  .  (/i;*kz),/z,  (G  .  R  ). 


(6.18a) 

(6.18b) 

(6.18c) 

(6.18d) 


(6.19a) 

(6.19b) 


Note  the  suppressed  dependence  of  cn  etc.  on  k  in  the  notation.  The 
remaining  terms  Efn  and  Df  are  zero  for  all  but  the  displacement  field, 


Eun  3  Evn  3  Ewn  3  0 

(6.20a) 

Du  =  Dv  =  Dw  =0 

(6.21a) 

Wz-zs>  ■-  rn(|kyl>ZlZs)  cn(ly>  1  co 

(6.20b) 

(ky,z,z$)  2  -  I'0(|ky|,z,zs)  /  cQ 

(6.21b) 

The  coefficients  (6.16)  through  (6.17)  have  all  been  constructed  so  as  to  be 
real.  The  above  Equations  (6.14)  through  (6.21)  represent  the  main  working 
results  for  numerical  computation  of  the  fields.  The  next  section  will 
consider  the  details  of  numerical  evaluation  of  the  various  terms  in  these 
equations. 
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The  above  results  hold  for  an  arbitrary  source  distribution  $0(x) 
subject  to  condition  (6.11).  If  the  source  has  symmetry  properties,  the 
above  results  simplify  considerably.  In  particular,  let  it  be  assumed  that 
the  source  distribution  is  symmetric  in  y  and  antisymmetric  in  x,  i.e. 

S0(-x,y)  =  -S0(x,y),  $0(x,-y)  =  S0(x,y). 

This  property  transfers  to  the  Fourier  transform  §0(kx,ky),  i.e.  it  also  is 
antisymmetric  with  respect  to  its  first  argument  and  symmetric  with  respect 
to  its  second.  It  then  follows  that  if  a  and  ky  are  both  real,  then  so  are 
both-  of  S0(ia,ky)  and  -1S0(a,ky).  Under  this  simplifying  assumption,  the 
above  solutions  may  be  put  in  the  common  form  : 


f  -  Af0(x) 


„S„  Mfn  Chn  e*P(-lx«nD 
n*o 


00 


*  „S,  !fn  Cin 


+  H(-x)  £ 
n=l 


2  Rfn  Crn  Afn(x)> 


(6.22) 


In  this  expression,  explicit  dependence  on  ky,  z,  z8,  etc.,  has  been 
suppressed  in  an  obvious  manner.  The  expression  is  valid  if  x>xm<lx  or  if 
x< xm, n *  The  first  pair  of  terms  are  present  both  ahead  of  and  behind  the 
source,  and  represent  a  localized  or  near-field  disturbance.  The  leading  A 
functions  are  defined  by 


Af0(x)  *  {-1,  -sgn(x) ,  -1,  -isgn(x)  },  f={f,w,u,v}.  (6.23) 

The  third  term  is  present  only  downstream  of  the  source  and  represents  the 
internal  wave  wake  generated  by  the  source.  Functions  Afn  are  defined  by 

Afn(x)  =  {  sin,  cos,  -sin,  icos}[x/in(ky)],  f={f,w,u,v}.  (6.24) 
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The  terms  C  are  related  to  the  source  function  by  : 

chn  *  S0[ii„(ky),k,].  C,„  »  S0[ii/„(ky),ky],  and 

5  •'  (6.25) 

and  are  real  on  account  of  the  above  note.  Thus  the  above  expression(6.22) 
is  real  and  even  in  ky  for  the  field  components  u,w,£,  and  is  imaginary  and 
odd  in  ky  for  field  v. 
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7.  NUMERICAL  CONSIDERATIONS 


(a)  Preliminaries 

Computation  of  a  physical  field  component  is  based  on  Equations 
(6.14)  through  (6.21).  The  integration  over  ky  is  implemented  efficiently 
using  well-known  Fast  Fourier  Transform  (FFT)  methods.  This  is  a  major 
reason  for  evaluation  of  the  kx-integrals  by  complex  variable  methods, 
leaving  ky  as  a  real  parameter,  and  it  essentially  allows  simultaneous 
computation  of  many  cross-track  points  at  each  fixed  x  and  z. 

The  heart  of  the  problem  is  therefore  the  computation  of  the 
eigenvalues,  eigenfunctions  and  various  coefficients  in  order  to  do  the 
modal  summations  of  (6.15b)  at  any  fixed  set  of  values  {x,ky,z,z,,c0}  with  a 
given  N2  profile. 

(b)  Differential  Equation  Solution;  General 

The  next  step  is  computation  of  values  of  the  Wronskian  A  in 
order  to  find  zeros.  This  entails  solution  of  the  homogeneous  differential 
equation  (3.4a),  together  with  boundary  conditions  (3.8),  for  a  given 
profile.  With  dependence  on  the  parameters  suppressed  for  now,  the 
differential  equation  can  be  written  as  a  first-order  linear  system, 

f(z)  «  A(z)  *(z),  (7.1a) 

where 

*  =  ^  =  Lq(z)  o  )  (7-lb) 

(  Note  the  convention  that  uppercase  bold  Greek  symbols  denote  2x1  column 
vectors,  and  uppercase  bold  underlined  Roman  symbols  denote  2x2  matrices.) 
The  solution  to  (7.1a)  can  always  be  expressed  in  the  form 


$(z)  =  F(z,z0)  i(z0) 


(7.2) 
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Here,  F  is  a  fundamental  matrix18  ("propagator")  relative  to  z0,  which 
satisfies  initial  conditions  F(z0,z0)=I2  (the  2x2  identity  matrix  )  and 
whose  columns  each  satisfy  (7.1a).  Note  the  basic  results 

F(z,z0)  =  F_1(z0,z),  (7.3a) 

£{ z.z,)  =  F(z,z#)  F(z0,z1).  (7.3b) 

The  determinant  of  F  is  independent  of  z  since  A  has  zero  trace,  and  the 
Initial  conditions  show  that 

det(F)  =  1.  (7.3c) 

It  follows  from  Equations  (3. 8), (3. 10),  (7.1)  and  (7.2)  that  the  Wronskian  A 
has  the  expression 

a  ~  E?? (z»0)F1?(z,d) ,  (7.4a) 

where  the  subscripts  refer  to  matrix  components.  The  left-hand  side  is 
independent  of  z,  and  so  the  right-hand  side  can  be  evaluated  at  any 
convenient  reference  depth  OizSd. 

If  it  should  prove  necessary  to  evaluate  derivatives  of  the 
Wronskian  with  respect  to  some  parameter,  two  approaches  are  possible.  If 
the  derivatives  of  F  are  available  analytically,  then  (7.4a)  rnay  be 
differentiated  directly  : 

A  »  *12(z,0)  F22(z,d)  ♦  F22(z,G)  *12(z,d) 

-  ^22(z,0)  F^Cz.d)  -  F22(z,0)  f12(z.d).  (7.4b) 

The  other  approach  is  to  differentiate  (3.10)  directly,  and  obtain  the 
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subject  to  the  appropriate  boundary  conditions  derived  from  (3.8). 


:c)  Layered  Profile 


Numerical  values  for  the  profile  Nz  will  generally  be  obtained 
from  measured  data  and  defined  on  a  discrete  mesh  of  points.  Many  methods 
can  be  employed  for  solution  of  the  differential  equation. 


One  possible  approach  is  a  three-point  finite  difference 

18  19  20 

scheme  *  ’  .  This  method  leads  to  a  tridiagonal  matrix  eigenvalue  problem, 


which  has  the  advantage  of  approximating  the  eigenvalues  of  the  differential 

4  ah  «  a  1  L  -  «  -  -  i  _  th  *  „  J  ,k  „  _  -i.  -  ♦  «  1  i  L  .  . .  _  L  1  I_  »  »  «  i 


equation  as  those  of  an  n  -order  matrix.  Although  the  matrix  elements  are 
extremely  simple,  an  excessively  fine  mesh  may  be  needed  to  adequately 
resolve  portions  of  the  eigenfunctions  with  large  vertical  wavenumbers. 


Another  possible  approach  is  to  use  piecewise  approximations  to 
the  profile  which  allow  the  differential  equation  to  be  solved  analytically. 
For  example,  piecewise  linear  approximations  lead21  to  Airy  function 
solutions  which  may  be  used  to  construct  the  fundamental  matrices  F  in  the 
various  segments.  While  such  methods  provide  accurate  solutions  on  a 
realistic  approximation  to  the  physical  profile  there  may  be  the  added 
burden  of  computation  of  the  various  higher  transcendental  functions  which 
arise. 


The  approach  to  be  used  here  is  to  use  the  simplest  possible 
approximation  to  the  profile,  namely  "layers"  within  which  N2  is  piecewise 
constant.  This  so-called  multilayer17,22  or  Thomson-Haskel l23  approach  leads 
to  eigenfunctions  which  are  piecewise  trigonometric  or  hyperbolic  elementary 
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functions.  This  approach  is  suitable24  for  the  present  problem  since  the 
basic  profile  data  is  N2  and  this  does  not  appear  in  differentiated  form  in 
the  differential  equation.  With  proper  treatment,  these  solutions  can  be 
suitably  stabilized  as  explained  below.  Note  that  this  approximation  may  be 
more  properly  thought  of  as  suitably  joined  local  WKB  approximations  to  the 
eigenfunctions  for  the  true  profile,  rather  than  a  stepped  approximation  to 
the  true  profile.  The  phase  of  the  local  trigonometric  solutions  allows  more 
efficient  eigenfunction  approximation  in  regions  of  large  vertical 
wavenumber  than  finite  difference  methods. 

Thus  for  a  given  profile,  a  set  of  points  {zo=0<z1<. . .<zn=d}  is 
introduced,  giving  n  layers,  of  which  the  mth  has  constant  value  Nm  and 
thickness  h^z^-z^.  An  example  of  a  such  profile  is  shown  in  Figure  2.  It 
is  based  on  measured12  summertime  temperature/salinity  data  taken  at  Ocean 
Station  P.  From  this  data,  N2-values  were  computed,  and  31  layers  fitted  to 
a  smoothed  version.  The  plot  displays  N(z)  in  units  of  cycles  per  hour  as  a 
function  of  depth  in  metres.  The  profile  shows  a  pronounced  pair  of  maxima. 
Only  the  top  1250m  of  data  is  plotted;  the  actual  profile  extends  to  a  depth 
d=4133m,  with  N=0  below  the  plotted  values. 

The  solution  form  (7.2)  for  any  points  z  and  z0  within  the  same 
layer  may  then  be  written  as 

*(Z)  ■  ICQ^iz-z.]  «U„).  (7.6) 

where  the  fundamental  or  propagator  matrix  T  is  defined  by 

cos(ah)  sin(ah)/a 

-a  sin(ah)  cos(ah) 

If  this  propagator  is  for  a  whole  layer,  it  is  denoted  T^,  i.e. 

X. 

Note  that  a  may  be  imaginary. 


(7.7) 


T(a,h)  ■ 


(7.8) 


Stability  Profile 


Figure  2.  The  layered  Brunt-Vaisala  frequency  profile  used  for  the  examples 
in  this  report.  It  is  based  on  a  smoothed  version  of  measured  data  taken  at 
Ocean  Station  P  during  the  summertime. 
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In  this  multilayer  Boussinesq  approximation,  appropriate 

conditions  on  the  solution  are  continuity  of  <f>  and  <f>‘  for  0<z<d  (this  is 
because  the  differential  equation  contains  no  derivatives  of  N2  ).  Thus  by 
application  of  the  general  form  (7.3b)  or  simple  iteration  of  (7.6),  it  is 
readily  seen  that  if  zr  is  within  layer  V,  then 

♦.(*.)  -ICw^r-Wo).  0-9«) 

and 

♦dUr)  ■  (7-9a) 

The  left-hand  sides  contain  the  information  necessary  for  evaluation  of  the 
Wronskian  A  at  an  intermediate  depth  zr.  In  principal,  this  reference  depth 
may  be  chosen  at  either  boundary,  and  so  the  Wronskian  can  be  expressed  as 

A  =  P12,  (7.10a) 

where  the  matrix  P  represents  the  fundamental  matrix  relative  to  z=0, 
evaluated  at  z»d,  i.e. 

£  3  Inln-1  *  “I,  •  (7-lOb) 

This  is  in  fact  the  first  equality  of  (3.10a)  in  the  present  notation. 

(d)  Auxiliary  Variables 

For  numerical  purposes  it  proves  convenient  to  work  in  terms  of 
new  parameters  (p,s)  in  place  of  (k2,c2).  These  are  defined  by 

p  »  5  *  Oc2  (7.11) 


A  convenient  combination  is  the  special  value 

Sc  3  ^raox7co» 


(7.12) 
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which  is  related  to  the  Froude  number4,  but  has  the  dimensions  of  reciprocal 
length  squared.  The  equations  inverse  to  (7.11)  express  (k,c)  in  terms  of 

(p.s), 

k2  >  sp,  c2  .  N2„x/S  (7.13) 

The  coefficient  Q  of  the  differential  equation  (3.4)  then  becomes 

Q(p,s;z)  =  s[v(z)-p],  (7.14) 

where  visa  dimensionless  profile 

K(z)  ■  N2(z)/nL.  (OSi/Sl).  (7.15) 

The  physical  eigenvalues  of  the  modified  form  of  (3.4)  then  define  a  set  of 
monotonical ly  increasing  curves  sm(p)  in  the  (p,s)  plane,  confined  to 
-“><p<l,  s^Q. 

Portions  of  the  first  21  oigencurves  sn(p)  in  Quadrant  I  of  the 
(p.s)-plane  are  shown  in  Figure  3.  These  correspond  to  the  profile  shown  in 
Figure  2,  and  are  displayed  in  a  log-log  (base  10)  format  in  order  to 
separate  the  various  curves  better.  Note  that  numbers  on  the  axes  give  the 
actual  values  of  p  ai.J  s,  not  their  logarithms. 

The  main  interest  is  in  those  points  which  lie  on  these 
eigencurves  and  also  satisfy  either  of  the  constraints  <u2  =>  c2  k2  in  the 
case  of  coefficients  (6.1«),  or  k2  »  0  in  the  case  of  (6.20b).  In  terms  of 
the  new  parameters,  the  first  constraint  translates  to 

(s-sc)p  »  k2,  (7.16a) 

and  the  second  to 

sp  =  k2.  (7.16b) 


t'Cy.  <jden-~  -*yC  !!'■?» fi&jfcv-.  ?■ - 


Figure  3.  Portions  of  the  first  21  eigencurves  in  Quadrant  I  of  the  (p,s) 
parameter  plane  for  the  profile  of  Figure  2.  Although  the  various  curves 
approach  each  other  with  exponentially  small  separation  at  a  discrete  number 
of  points,  they  do  not  intersect  each  other.  The  illusion  of  two  sets  of 
curves  is  due  to  the  presence  of  the  two  peaks  in  the  profile,  as  discussed 
in  the  text.  Also  shown  are  three  constraint  curves  along  which  eigenvalues 
are  sought  for  a  fixed  value  of  kv  and  source  speed. 
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Both  represent  hyperbolae  in  the  (p,s)-plane,  and  have  the  s-axis  as  a 
vertical  asymptote.  The  second  constraint  hyperbola  intersects  the 
eigenvalue  curves  only  in  the  first  quadrant  p,s^Q,  providing  a  denumerable 
infinity  of  intersections  {pm* s,,,}  corresponding  to  real  eigenvalues.  The 
first  constraint  has  a  horizontal  asymptote  in  s>0,  and  so  provides  a  first 
denumerable  set  of  solutions  {p”,s“}  in  the  second  quadrant  s*0,p£0 
corresponding  to  the  imaginary  eigenvalues  of  (6.19a),  and  a  second  set 
{p*,s*}  in  the  first  quadrant  s^O.p^O  corresponding  to  the  real  eigenvalues 
of  (6.19b).  The  effect  of  a  increase  in  source  speed  c0  is  to  shift  the 
horizontal  asymptote  in  (7.16a)  downwards,  while  the  effect  of  an  increase 
in  ky  is  to  increase  the  distance  of  the  hyperbolae  from  the  origin.  The 
eigencurves  themselves  are  a  property  purely  of  the  N2-profile.  The 
imposition  of  the  c0  and  ky  values  merely  dictates  how  these  curves  are 
"sampled"  in  the  integrations  leading  to  the  physical  fields. 

For  each  fixed  ky  (to  be  dictated  by  the  integration  scheme 
chosen  for  the  ultimate  integration  of  (6.14a)  ),  the  finding  of  eigenvalues 
is  thus  conveniently  reduced  to  a  one-dimensional  problem  of  searching  along 
the  hyperbolae  (7.16),  which  can  be  expressed  in  terms  of  a  single  suitable 
variable.  Hyperbola  (7.16a),  for  example,  may  be  described  by  a  parameter 
'V  by  way  of  the  definitions 

p(v)  »  o(v)-v/2;  s(v)  =  s0+o(v)+v/2  (Quadrant  I),  (7.17a) 

p(v)  »  -<r(v)-v/2;  s(v)  *  s0-«r(v)+v/2  (Quadrant  II),  (7.17b) 

where 

cr(v)  a  (  kj  ♦  v2/4  )1/2,  and  -«<v<®;  (7.18) 

the  signs  in  (7.17)  have  been  chosen  so  as  to  make  s * ( v) >0 ,  p‘(v)<0  on  both 
arms  of  the  hyperbola.  Hyperbola  (7.16b)  may  be  parametrized  by  equations 
analogous  to  (7.17)  with  sc  replaced  by  0. 

Sample  hyperbolae  of  the  form  (7.17)  for  several  ky-values  are 
also  shown  on  Figure  3.  In  the  log-log  form  used,  these  hyperbolae  become 
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almost  linear;  log(p)  +  1og(s-sc)  =  21og(ky).  The  speed  used  for  these 
constraint  illustrations  was  co=0.5  ms-1,  and  the  ky-va1ues  correspond  to 
{ j (rr/3072)  rad  nf1;  j=  25,  271,  513  }.  The  intersections  of  any  one  of  these 
hyperbolae  with  the  eigencurves  give  the  corresponding  set  of  modal  points 
{ P» » s«}  described  above,  for  the  value  of  ky  under  consideration. 

Thus  once  c0  and  ky  are  fixed,  the  constrained  eigenvalue  problem 
has  only  the  single  parameter  v,  and  involves  finding  zeros  of  the  (real) 
function 

D(v)  ■  A[p(v),s(v)]  (-»<v<»),  (7.19) 

where  A  is  the  Wronskian  (3.9),  expressed  in  terms  of  the  new  parameters. 
Figure  4  shows  plots  of  D(v)  along  the  three  constraint  hyperbolae  plotted 
of  Figure  3.  Note  that  the  v-ranges  on  each  plot  are  different;  the 
parametrization  (7.17)  varies  with  ky. 

(e)  Effects  of  Interior  Evanescent  Regions 

In  practice,  the  presence  of  evanescent  layers  (  imaginary 
vertical  wavenumbers  Qa<0  ),  for  certain  values  of  the  parameters  (p,s),  is 
a  well-known  source  of  difficulty11,25,28.  As  may  be  seen  from  (7.14),  the 
differential  equation  coefficient  for  the  mth  layer  is  Qraas(i/m-p) .  As  noted 
above,  the  imaginary  eigenvalues  have  -®<p£.0,  s>0  and  so  all  layers  are 
oscillatory.  However  the  real  eigenvalues  have  s>0  and  0£p£l.  Thus  for  a 
given  value  of  p,  all  layers  with  i/m<p  will  be  evanescent.  If  not  properly 
handled,  these  regions  may  cause  numerical  exponential  overflows  and 
inappropriate  approximations  for  eigenfunctions.  Moreover,  such  regions  may 
cause  eigenvalues  to  be  exponentially  close  together  and  so  cause  multiple 
eigenvalues  in  a  numerical  implementation.  Such  multiplicities  and  overflows 
must  be  handled  carefully,  but  this  can  be  done  in  a  straightforward  manner 
as  will  be  explained  later.  Before  a  general  discussion  however,  it  is 
instructive  to  consider  a  symmetric  three-layer  profile. 
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Figure  4.  Plots  of  the  scaled  Wronsk’ 
on  Figure  3.  The  three  values  of  ky  cf 
the  close  separation  of  the  fourth 
causes  their  apparent  coalescence  intc 


)  along  the  three  cuts  depicted 
lcrease  down  the  page.  Note  that 
fth  zeros  in  the  central  plot 
ole  zero. 
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(f)  A  Symmetric  Three-layer  Example 

The  symmetric  three-layer  profile  is  shown  in  Figure  5. 


z 

Figure  5.  Symmetric  three-layer  profile.  Note  that  the  origin  has 
been  shifted  for  purposes  of  this  discussion. 

It  is  convenient  for  purposes  of  discussion  to  consider  the  surface  to  be  at 
z=-a-b  and  the  bottom  at  z=a+b.  Thus  the  centre  layer  has  thickness  '2b1, 
and  the  differential  equation  coefficient  there  is  /S2  s  s(i/-p).  The  two 
flanking  layers  each  have  thickness  'a'  and  coefficient  ct2  =  s(l-p).  Within 
the  centre  region  the  general  solution  necessarily  has  the  form 

$(z)  =  A  cos(/?z)  +  B  sin(0z)/0. 

It  may  then  be  shown  that 

0(±a±b)  =  LA  ±  MB, 


where 


L  s  cos(/?b)cos(aa)  -  (0/a)sin(/Sb)sin(aa) , 


M  e  sin(/3b)cos(aa)//?  +  cos(/Sb)sin(aa)/a. 
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The  boundary  conditions  0(± a±b)  =0  then  result  in  the  homogeneous  equations: 

LA+MB=G,  LA-MB=Q, 
so  that  either 

1=0,  B=0  (even  eigenfunctions) , 
or 

M=0,  A=0  (  odd  eigenfunctions). 

In  terms  of  the  standard  notation,  it  may  be  shown  that  for  this  problem 
A  =  2LM, 

i.e.  the  Wronskian  factors  into  a  product  of  two  terms.  The  eigenvalue  , 
conditions  are  thus 

acot(aa)  =  0tan(/9b)  (even), 
or 

acot(aa)  =-0cot(jSb)  (odd). 

If  the  eigenfunctions  are  normalized  to  unit  slope  at  the  upper  (z=-a-b) 
boundary,  so  that  they  correspond  to  the  "upper"  solutions  <f>0  of  the 
standard  notation,  they  may  be  written  in  the  form: 


-b-a<z<-b, 

-b<z<b, 

b<z<a+b. 


where  the  subscripts  here  denote  five n  or  Odd  eigenfunctions  respectively. 

If  i/<p<l,  then  the  centre  layer  is  evanescent,  and  <9=1101.  The  outer  layers 


-47- 


are  oscillatory,  and  the  Wronskian  A  may  be  written  as: 

A  »  [$inh(2l/?lb)sin(aa)/(al/?! )]  x 

[acot(aa)  +  ljS|coth(  1  ^ i  b]  [acot(aa)  +  !jSltanh(  l)S|b]. 

The  above  eigenvalue  conditions  reduce  to 

acot(aa)  »-l0ltanh( l/Sib)  (even), 
or 

acot(aa)  a-|/?lcoth(l/J|b)  (odd). 

If  the  product  l/?!b  is  large,  two  points  may  be  noted.  Firstly,  the  leading 
sinh(2l/SIb)  causes  the  Wronskian  A  to  become  exponentially  large  except  near 
its  zeros.  Secondly,  both  the  t anh ( l/Slb)  and  coth  (l/Jlb)  right  hand  sides 
approach  the  common  value  1,  and  so  the  eigenvalues  merge  into  pairs  with 
exponentially  small  separation.  These  factors  exhibit  both  the  exponential 
overflow  and  multiple  zero  problems  mentioned  in  the  preceding  section.  As 
will  be  shown  below,  the  exponential  amplitude  factor  is  easily  handled,  but 
the  presence  (in  a  numerical  sense)  of  multiple  zeros  may  cause  problems. 

Consider  now  the  eigenfunctions  in  the  presence  of  this 
evanescent  layer.  In  this  centre  layer  the  even  eigenfunctions  have  the  form 

0E(z)  *  o_1sin(aa)cosh(l|?lz)/cosh(ljSlb), 

while  the  odd  eigenfunctions  have  the  form 

4>0{  z)  =  -a_1sin(aa)sinh(  l(?lz)/sinh(  l/Slb) . 

For  large  IjSlb  products,  the  eigenfunctions  in  the  outer  layers 
corresponding  to  two  closely-separated  modes  may  be  numerically 
Indistinguishable,  apart  from* a  possible  change  in  sign  in  one  of  the 
oscillatory  layers.  This  may  be  seen  from  the  above  general  expressions  for 
the  eigenfunctions,  where  an  even  and  adjacent  odd  eigenfunction  may  have 
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values  of  a  which  are  exponentially  close  together.  In  the  evanescent  layer, 
the  eigenfunctions  are  essentially  zero  over  most  of  the  layer.  The 
additional  interior  zero  of  the  sinh  term  ensures  that  in  principal  it 
represents  the  next  higher  mode,  but  in  practice  it  may  be  impossible  to 
distinguish  numerically  between  the  exponentially  small  values  of  the  cosh 
term  and  the  "true"  zero  in  the  sinh  term. 

The  exponential  smallness  in  the  interior  evanescent  layer  at 
first  glance  appears  to  represent  “tunnelling"  of  a  finite  amplitude 
eigenfunction  through  a  region  of  extreme  smallness,  only  to  reappear  with 
finite  amplitude  in  the  other  layer.  The  resolution27  lies  in  consideration 
of  the  eigenfunction  expansion,  including  the  source  term.  Suppose  that  the 
source  is  in  one  of  the  outer  oscillatory  layers  (cf  sound  channels  in 
acoustics  ).  By  grouping  the  modal  summation  in  pairs  over  adjacent  modes  in 
those  terms  where  the  eigenvalues  coalesce,  it  can  be  seen  that  the  source 
excites  these  two  modes  equally  in  its  own  layer,  but  the  sign  change  of  the 
eigenfunctions  in  the  opposite  layer  removes  any  contribution  there.  Thus 
for  numerical  purposes,  a  modified  set  of  eigenfunctions  may  be  defined  by 
doubling  the  amplitude  in  one  half  of  the  profile  and  zeroing  it  in  the 
other.  Thus  when  an  eigenvalue  pair  is  numerically  indistinguishable,  It  is 
counted  twice,  and  used  to  form  two  numerically  orthogonal  eigenfunctions 
which  are  disjoint  and  normalized  over  half  the  profile  (which  effectively 
doubles  their  amplitude  in  each  layer  in  keeping  with  the  pairwise 
regrouping  in  the  eigenfunction  expansion).  It  is  then  the  source  which 
dictates  which  of  the  disjoint  (modified)  eigenfunctions  is  excited.  The  end 
result  is  that  energy  excited  in  a  given  "waveguide"  subsection  of  the 
overall  profile  remains  confined  to  that  region. 

(q)  Treatment  of  Interior  Evanescent  Regions 

Now  return  to  consideration  of  the  more  general  problem  in  which 
there  may  be  multiple  interior  evanescent  regions.  The  exponential  growth 
problem  can  be  easily  circumvented  simply  by  factoring  out  an  exponential 
term  from  the  matrix  for  each  evanescent  layer  in  the  product  (7.10). 
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Specifically,  (7.10)  may  be  recast  in  the  form 


where 


P  =  9  exp[  /J  ImVQ(z')  dz'  ], 


2  =  E^En-r-Ei. 


a,  s  X,  exp[-Im(haVQra)]. 


(7.20a) 


(7.20b) 


(7.20c) 


If  the  layer  is  oscillatory  (Q>0)  ,  then  the  matrix  R  is  identical  to  T  as 
defined  in  (7.7).  However,  if  the  layer  is  evanescent,  then  R  has  the  form 


[l+exp(-2ah)]/2  [l-exp(-2ah)]/(2o) 

a[l-exp(-2ah)]/2  [l+exp(-2ah)]/2 


(7.21) 


where  here  asabs(ImVQ).  The  elements  of  each  matrix  R  are  bounded,  and  for 
an  evanescent  layer  they  are  also  positive.  Moreover,  if  several  adjacent 
interior  layers  are  evanescent,  then  the  product  of  the  corresponding  R's  is 
also  bounded  and  has  positive  elements.  There  is  then  no  numerical 
difficulty  (or  "instability")  in  evaluation  of  the  matrix  product  (7.20b). 

nc  no 

The  above  approach  is  similar  to  ones  used  elsewhere  '  . 
Specifically,  (discontinuous)  amplitude  control  functions  x0.d  f°r  solutions 
of  the  differential  equation  are  defined  by  way  of 


X0(z)  a  exp[  f*  ImVQ(z')  dz'  ], 
Xd(z)  s  exp[  jt  ImVQ(z')  dz1  ]; 


(7.22a) 


(7.22b) 


the  integrands  here  are  by  convention  either  zero  or  positive.  Scaled 
“eigenfunctions"  t0  (S  are  then  defined  by 


*0(z)  a  X0(z)lo(z)'’  §d(z)  3  Xd(z)la(z)- 


(7.23a,b) 


The  functions  x  limit  any  exponential  growth  in  the  scaled  eigenfunctions 
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away  from  their  respective  boundary  and  the  modified  functions  £  are 
bounded. 

On  account  of  (7.20a),  a  zero  of  §  is  one  of  P,  and  hence  an 
eigenvalue.  The  eigenvalue  search  is  then  reduced  to  finding  zeros  of  the 
(1,2)  element  of  the  bounded  matrix  Q.  It  should  be  noted  that  extensive 
numerical  experimentation  verifies  that  there  is  no  need  to  evaluate  the 
Wronskian  at  an  interior  reference  depth  for  the  purposes  of  root-finding. 
However,  as  will  be  discussed  later,  more  care  is  needed  with  construction 
of  the  eigenfunctions,  where  the  evanescent  layers  may  introduce  artificial 
amplification  of  the  eigenfunctions  in  inappropriate  regions  of  the  profile. 

It  may  be  noted  that  for  an  evanescent  layer,  det(R)  =  exp(-2ah). 
Thus  for  evanescent  layers  with  large  a-h  products,  or  for  a  set  of  adjacent 
layers  whose  integrated  a-h  product  is  large,  the  determinant  of  the 
corresponding  matrix  can  vanish  (underflow)  numerically.  This  feature  in 
effect  factors  the  eigenvalue  problem  into  two  subproblems  in  which  an 
artificial  rigid  boundary  (or  appropriate  radiation  condition)  is  introduced 
into  the  evanescent  region.  This  is  the  generalization  of  the  exact 
factorization  which  arises  in  the  symmetric  three-layer  model  above. 
Specifically,  suppose  that  the  profile  contains  an  interior  region  of  high 
evanescence.  Then  the  factors  of  (7.20b)  may  be  regrouped  as 

fi  *  LEU,  (7.24a) 

where  L  denotes  the  product  of  matrices  for  all  layers  from  below  the 
evanescent  region  to  the  bottom,  E  represents  all  layers  within  the 
evanescent  region  under  consideration,  and  U  represents  all  layers  above  the 
evanescent  region.  (Note  that  L,y  may  themselves  contain  evanescent  regions, 
in  which  case  the  present  argument  is  Iterated).  Numerically,  det ( E) ,  so 
i^22~^i 2^21  *  Substitution  of  this  relationship  into  (7.24a)  leads  to  the 
result 


fli2  “  (ill)  12(^)12^12* 


(7.24b) 
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The  denominator  is  strictly  positive.  The  two  numerator  factors  in  effect 
represent  separated  eigenvalue  problems.  In  the  first  sub-problem,  the 
“upper"  profile  (represented  by  U  )  overlies  the  evanescent  region  (E),  and 
In  the  second,  the  evanescent  region  overlies  the  "lower"  structure 
represented  by  L.  The  separated  problems  thus  correspond  to  the  total 
propagator  matrices  IE  and  EU.  Each  factor  will  contain  its  own  zeros, 
leading  to  the  possibility  of  multiple  zeros  in  the  left  hand  side  of 
(7.24b),  depending  on  the  details  of  the  profile.  In  general  these  will 
occur  only  at  discrete  points,  where  families  of  eigencurves  corresponding 
to  each  subsection  approach  each  other  exponentially  closely  and  in  effect 
appear  to  cross. 

This  effect  is  clearly  visible  in  Figure  3.  It  should  be 
emphasized  that  the  eigencurves  do  not  intersect  each  other  -  they  do 
however  approach  each  other  with  exponentially  small  separation  at  a 
discrete  number  of  points.  These  points  of  close  approach  give  the  illusion 
of  two  distinct  sets  of  curves  for  larger  (p,s)  values.  These  pseudo-curves 
in  fact  correspond  to  eigencurves  appropriate  to  the  two  regions  on  either 
side  of  the  interior  minimum  of  the  profile  of  Figure  2,  which  is  evanescent 
for  sufficiently  large  values  of  p. 

The  interference  effect  of  the  two  factors  in  (7.24b)  is  also 
evident  in  Figure  4.  It  is  more  pronounced  for  higher  values  of  ky,  when  the 
interior  evanescence  is  higher.  The  central  plot  of  Figure  4  exhibits  two 
very  closely  separated  zeros  corresponding  to  modes  4  and  5.  This  is  also 
evident  in  Figure  3  where  the  corresponding  constraint  curve  appears  to  pass 
through  an  intersection  of  the  fourth  and  fifth  eigencurves. 

The  presence  of  highly  evanescent  Interior  regions  thus  need  not 
cause  problems  In  searching  for  eigenvalues,  but  a  complete  numerical  scheme 
must  be  able  to  find  and  properly  interpret  multiple  roots.  Note  that  a  more 
sophisticated  scheme  could  in  effect  avoid  multiple  root-finding  by 
appropriately  subsectioning  the  profile,  depending  on  the  range  of 
parameters,  and  looking  for  separate  families  of  eigenvalues  appropriate  to 
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each  section,  as  defined  by  (7.24b).  These  eigenvalues  could  then  be 
appropriately  sorted  later.  This  method  would  correspond  to  finding  zeros  in 
the  separated  factors  as  in  (7.24b)  when  evanescence  is  significant,  but 
involves  additional  'bookkeeping'  for  the  various  regions.  Note  that  in 
terms  of  the  parameters  (p,s),  the  physical  eigenvalues  all  have  s>0,  so 
that  regions  of  evanescence  are  controlled  only  by  the  parameter  p.  A 
subsectioning  approach  would  then  be  easiest  to  implement  for  problems 
requiring  eigenvalues  at  fixed  values  of  p  (i.e.  <e2  ).  For  the  moving  source 
problem  however,  both  p  and  s  are  constrained  to  a  hyperbola,  and  so  the 
number  and  extent  of  evanescent  regions  varies  with  the  parameter  v,  thereby 
presenting  additional  complications  in  implementing  a  subsectioning  scheme 
of  this  type.  Such  an  approach  was  not  found  necessary  for  finding 
eigenvalues  in  the  profiles  used  to  date,  however. 

(h)  Computation  of  the  Eigenfunctions 

As  noted  above,  interior  evanescent  layers  may  induce  numerical 
multiple  zeros  for  certain  values  of  the  parameters  (p,s).  The  problem,  and 
common  cause  of  "instability"  in  the  presence  of  evanescent  regions,  comes 
in  computation  of  the  corresponding  eigenfunctions. 

For  example,  in  the  three  layer  profile  above,  it  will  prove 
difficult  to  evaluate  the  "true"  eigenfunctions  by  integration  in  one 
direction  due  to  the  exponential  smallness  in  the  centre  region.  It  Is 
easier  to  compi’U  the  "modified"  eigenfunctions  described  above  however,  by 
integrating  Inward  from  each  boundary  to  a  reference  depth  in  one 
oscillatory  layer,  and  then  repeating  the  procedure  for  a  second  reference 
depth  in  the  other  oscillatory  layer. 


Thus  in  the  general  profile,  the  presence  of  Interior  evanescent 
regions  may  cause  artificial  amplification  in  construction  of 
eigenfunctions,  if  the  reference  depth  is  not  chosen  appropriately. 
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A  "standard"7  N2  profile  has  a  single  peak  (thermoel ine)  at  some 
depth.  The  N2  tapers  to  zero  toward  the  mixed  surface  layers  and  also  toward 
the  bottom.  There  are  no  multiple  eigenvalues  for  such  a  profile,  and 
computation  of  the  eigenfunctions  is  easily  handled25  by  choosing  a 
reference  depth  zr  in  the  layer  containing  the  maximum  N2-value.  The  scaled 
eigenfunctions  f  are  evaluated  from  the  surface  down  to  zr  and  from  the 
bottom  up  to  zr.  These  are  then  suitably  matched  at  z=zr,  and  finally,  the 
exponential  factors  (7.22)  are  included  to  produce  exponential  decay  through 
any  evanescent  layers  toward  the  boundaries. 

In  the  presence  of  interior  evanescent  regions,  the  above 
procedure  must  be  modified,  and  the  definition  of  the  eigenfunctions 
modified  as  was  done  in  the  case  of  the  symmetric  three-layer  profile  above. 
In  parameter  ranges  where  the  zeros  of  the  Wronskian  are  well-separated, 
there  is  no  problem.  The  matrix  product  (7.20b)  for  the  current  eigenvalue 
can  be  split  into  subproduets  over  regions  separated  by  evanescent  layers 
according  to  some  prechosen  cutoff.  (This  is  easily  accomplished  by 
examining  Ira{Qn}  for  each  layer.  )  Specifically,  for  each  group  of  adjacent 
interior  evanescent  layers,  a  depth  zp  may  be  defined  at  the  "centre  of 
evanescence".  Successive  such  points  then  decompose  the  profile  into 
alternating  regions,  each  comprising  a  central  oscillatory  section,  bounded 
on  each  side  by  an  evanescent  one.  Thus  the  product  Q  may  be  symbolically 
written  in  the  form 

9  23  Bi  Si  •••  Em.  (7.25a) 
with  each  term  having  the  symbolic  form 

R'  «  EOE.  (7.25b) 

where  the  subscripts  E  and  0  here  indicate  evanescent  and  oscillatory;  one 
of  the  terms  E  may  in  fact  be  the  identity  matrix  near  a  boundary.  One  or 
more  of  these  matrices  will  have  a  zero  in  its  (1,2)  entry,  if  the 
partitioning  cutoff  is  numerically  large  enough  to  induce  factoring  as  in 
(7.24b). 
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If  the  zero  is  a  simple  one,  then  a  reference  depth  within  the 
(unique)  subsection  with  an  (approximate)  zero  in  its  (1,2)  entry  is  chosen. 
Construction  of  the  eigenfunctions  then  proceeds  as  in  the  case  of  a  single 
thermocline  by  integration  of  the  scaled  eigenfunctions  f  to  this  chosen 
reference  depth.  After  matching  at  the  reference  depth,  inclusion  of  the 
exponential  "tails"  represented  by  the  functions  x  over  the  intervals  (0,zr) 
and  (zr,d)  effectively  zeros  contributions  in  other  regions. 

In  the  case  of  an  eigenvalue  of  multiplicity  'm',  there  will  be 
'm'  such  regions.  The  eigenvalue  is  counted  m  times,  and  m  (disjoint) 
eigenfunctions  constructed  from  it  by  repeating  the  above  procedure  for  a 
reference  depth  chosen  once  within  each  of  the  m  regions  containing  zeros. 

The  net  effect  is  to  produce  an  eigenfunction  which  pertains  to 
the  subsection  in  question  and  has  exponentially  small  (numerically  zero) 
amplitude  elsewhere.  The  procedure  thus  automatically  induces  a  modified 
eigenfunction  expansion  as  discussed  for  the  example  above,  and  properly 
accounts  for  channelling  of  energy  in  interior  ducts;  it  is  the  source  term 
which  selects  the  appropriate  eigenfunction  combination. 

Figure  6a  shows  the  first,  eleventh  and  twenty-first 
eigenfunctions  over  the  upper  4000m  of  the  profile  of  Figure  2.  These  have 
been  normalized  so  that  =  N^ax  (see  Equations  (3.6a)  and  (7.30a)  ).  The 
corresponding  values  of  ky  and  c0  are  (rr/768)  rad  nf1  and  0.5  ms"1 
respectively.  For  such  a  small  wave  number  there  are  no  significant 
evanescent  regions  except  for  the  bottom  section.  These  long  wavelength 
components  'see'  the  entire  profile. 

The  effects  of  interior  evanescence  and  eigenfunction 
modifications  are  considered  for  modes  4  and  5  near  the  central  ky  value 
illustrated  in  Figures  3  and  4.  This  central  kyj  =  jn/3072,  j=271, 
corresponds  to  a  pair  of  very  closely  separated  zeros  in  the  Wronskian  D(v). 
Although  it  is  not  apparent  from  the  Figures,  the  two  zeros  are  in  fact 
easily  resolved  numerically,  yet  graphically  Illustrate  the  eigenfunction 
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behaviour.  Figure  6b  shows  Mode  4  corresponding  to  j=269  and  j=273,  while 
Figure  6c  presents  Mode  5  for  the  same  pair  of  ky  values.  The  two 
eigenfunctions  exhibit  the  opposite  behaviour  in  the  upper  and  lower 
oscillatory  sections  of  the  profile. 

With  reference  to  Figure  3,  it  is  evident  that  the  point  of 
closest  approach  of  modes  4  and  5  for  the  full  profile  corresponds  to  the 
first  mode  of  the  lower  oscillatory  subsection  of  the  profile  (i.e.  the 
first  apparent  'curve'  with  vertical  asymptote  near  p*0.8),  and  to  the 
fourth  mode  of  the  upper  oscillatory  region.  Thus  the  eigenfunctions  would 
be  expected  to  have  three  zeros  tn  the  upper  section  and  no  zeros  in  the 
lower  one,  as  is  evident  in  Figures  6b  and  6c. 

In  Figure  6b,  the  upper  plot  is  for  a  ky-value  slightly  to  the 
left  of  the  'intersection'  point.  As  is  evident  from  Figure  3,  this 
corresponds  to  mode  1  of  the  lower  section,  and  so  is  computed  for  a 
reference  depth  there,  and  decays  exponentially  toward  the  boundaries.  The 
lower  plot  of  Figure  6b  is  for  a  slightly  higher  ky,  to  the  right  of  the 
'intersection'  point  and  so  strongly  excites  mode  4  in  the  upper  section, 
where  the  reference  depth  was  chosen. 

In  Figure  6c,  the  opposite  behaviour  is  exhibited  for  mode  5.  The 
only  difference  is  the  sign  change  in  the  lower  section,  which  is  due  to  an 
additional  zero  within  the  evanescent  region,  and  accounts  for  this 
representing  mode  5  for  the  full  profile.  However,  for  the  purposes  of 
evaluating  a  fluid  dynamical  field  in  which  the  modes  are  summed  and 
integrated  after  inclusion  of  a  source  term,  it  is  evident  that  the  'label' 
(mode  number)  applied  to  a  particular  mode  is  of  little  importance.  A  source 
in  the  lower  section,  for  ky  values  near  the  above  'intersection'  point, 
excites  mode  1  of  the  lower  section,  regardless  of  the  fact  that  this  is 
labelled  either  mode  4  or  mode  5  for  the  full  profile  for  adjacent  ky’s. 

For  purposes  of  illustration,  Figure  6d  shows  both  modes  4  and  6 
plotted  together  for  the  above  pair  of  ky  values.  Due  to  the  exponential 
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Figure  6b.  Normalized  eigenfunction  for  mode  4,  for  two  closely  separated 
intermediate  values  of  ky,  one  on  either  side  of  the  central  constraint 
curve  shown  on  Figure  3.  The  presence  of  interior  evanescence  causes  the 
sudden  behaviour  change.  For  the  upper  value,  the  eigenfunction  looks  like 
mode  1  appropriate  to  the  lower  section  of  the  profile  in  isolation.  For  the 
lower  plot,  it  looks  like  mode  4  appropriate  to  the  upper  section.  This 
behaviour  can  be  expected  from  Figure  3. 


Figure  6d.  Figures  6b  and  6c  plotted  together  to  illustrate  the  possibility 
of  ambiguity  or  instability  in  a  standard  shooting  method  for  eigenfunction 
computation.  Both  the  eigenfunction  and  its  derivative  are  exponentially 
small  in  the  interior  region,  and  curves  looking  like  any  one  of  the  6  plots 
of  Figures  6b  through  6d  may  be  acceptable  numerical  solutions  near  closely 
separated  eigenvalues. 

smallness  of  the  eigenfunctions  in  either  the  upper  or  lower  sections  of  the 
profile,  these  two  curves  are  approximately  what  is  obtained  from  the  sum 
and  difference  of  modes  4  and  5  near  the  intersection  point.  If  the 
evanescence  were  sufficiently  large  for  the  eigenvalue  to  be  a  double  one  in 


-60- 


a  numerical  sense,  then  a  pair  of  eigenfunctions  looking  like  the  upper  and 
lower  curves  of  Figure  6d  would  be  an  acceptable  solutions,  the  upper  one 
connected  by  a  cosh  term  and  the  lower  one  by  a  sinh  term  in  the  evanescent 
region,  precisely  as  occurred  for  the  symmetric  three-layer  example  above. 
For  numerical  purposes  however,  it  preferable  to  work  with  the  modified 
forms  of  Figures  6b  and  6c. 


(i)  Connection  with  Global  Matrix  Methods 

There  is  a  close  connection  between  the  "propagator"  approach  for 
solving  the  differential  equation  stepwise,  and  "global  matrix"  methods11. 
Specifically,  if  0,  denotes  the  amplitude  of  a  solution  of  the  differential 
equation  at  the  ith  interface  (i=Q...n),  then  it  may  be  shown  from  the 
propagator  equations  (by  eliminating  the  interface  derivative  values  <f>‘ ) 
that  the  interface  values  satisfy  the  three-term  recurrence  relation: 


where 


and 


ai  0i-i  +  (bi+bi+i)  0.i  +  ai+i  0i+i  s  0, 

a,  a  qj/sinCqjh,) , 

b,  s  +qicos(qihi)/sin(q,h]) , 


(7.26a) 

(7.26b) 

(7.26c) 

(7.26d) 


Together  with  the  boundary  conditions  <j>0  =  <f>n  =  0,  the  above  is  equivalent 
to  a  homogeneous  tridiagonal  (n-l)x(n-l)  matrix  system: 
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(7.27a) 


which  may  be  written  concisely  in  matrix  form  as 


C-y  -  0,  (7.27b) 

where  the  matrix/vector  notation  here  refers  to  dimension  n-1.  For  a  non¬ 
trivial  solution  det(C)=0.  In  the  presence  of  an  interior  highly  evanescent 
layer,  one  of  the  a's  in  (7.27a)  will  be  exponentially  small  and  so  may  be 
numerically  zero  in  a  computer  implementation.  In  this  case,  the  system 
(7.27b)  may  be  written  in  partitioned  matrix  form  as: 
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In  this  case,  there  is  again  a  factorization  in  the  eigenvalue 
since 


(7.28) 

criterion, 


det(C)  ■  det(Ct)  det(C2).  (7.29) 

The  coupling  between  the  pairs  in  the  partition  is  absent.  If  det(C1)=0, 
then  a  numerically  acceptable  solution  to  (7.29)  has  the  form  (y^O),  and 
similarly  if  det(C2)=0  a  solution  has  the  form  (0,y2).  In  the  case  of  a 
coincident  pair  of  eigenvalues,  the  interpretation  is  once  again  a  pair  of 
disjoint  eigenfunctions  {(y1»0),(0,y2)},  the  selection  of  which  is  again 
made  by  the  source  term.  An  alternative  pair  of  acceptable  solutions  is  any 
pair  of  linearly  independent  linear  combinations  of  the  above  pair,  and  in 
particular  the  sum  and  difference,  namely  ( (yt ,±y2) } . 


■m 
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(i)  Integrands  for  the  Fields 

The  preceding  discussion  has  explained  how  eigenpairs  (pra,sra)  may 
be  computed  for  a  specified  set  of  ky  values,  and  how  a  corresponding  set  of 
(possibly  unnormalized)  eigenfunctions  may  be  computed  in  the  modified  sense 
described  above.  The  next  step  is  construction  of  the  various  terms  needed 
for  evaluation  of  the  fields,  as  given  by  Equations  (6.14a)  and  (6.15) 
through  (6.21). 

The  normalization  factors  (3.6a)  are  most  easily  computed 
directly  from  the  definition  by  integration: 

«  ■  C,  St  dz.  (7.30a) 

The  integration  here  is  done  by  breaking  it  into  a  sum  over  individual 
layers  and  integrating  the  various  terms  (which  are  simple  trigonometric  or 
hyperbolic  functions)  analytically.  It  is  also  easy  to  evaluate  the  related 
integrals 


4  -  It  £(I)  dz  (7.30b) 

simultaneously.  The  slopes  s^(p)>0  of  the  eigencurves  may  be  expressed  in 
terms  of  these  integrals.  Repetition  of  the  calculations  which  lead  to 
(6.4),  in  terms  of  the  new  parameters  (p,s),  leads  to  the  equation 

ft  C  sm(P)  (‘'-P)  -  PSra(p)  ]  <t>l  dz  =  0, 
from  which  it  is  seen  that 

-  P  *1  ]•  (7.31) 

The  phase  velocities  and  horizontal  total  wavenumbers  follow  from  (7.13) 

k2  =  sp,  c2  .  N2ox  /  s  (7.32) 
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Values  of  kx  are  either  zero  or  are  obtained  from 


<1  ■  psc.  (7-33) 

depending  on  the  field  integrand  under  consideration.  Values  of  kx  and  k 
will  be  imaginary  for  p<0.  The  other  quantity  of  interest  is  the  derivative 
(6.4),  which  in  the  present  variables  translates  to 

daj/dk2  *  NjJax  /  [s.(p)  +  p  s;(p)  ].  (7.34) 

With  these  results,  and  a  description  of  the  source  transform  $0,  all  of  the 
terms  in  the  integrands  (6.15)  may  be  computed. 

Figure  7  shows  the  first  21  real  eigenvalues  kxni(ky), 
corresponding  to  (7.33).  These  correspond  to  the  singular  curves  in  the 
(kx,ky)  plane,  which  were  of  concern  in  evaluation  of  the  inverse  Fourier 
transform  (5.4).  These  curves  are  again  for  the  profile  of  Figure  2,  but 
with  a  source  speed  of  2.5ms"1.  For  larger  k-values,  there  are  again  clearly 
two  families  of  curves,  which  give  the  illusion  of  intersection, 
correspond!  j  to  the  two  peaks  in  the  profile.  For  each  N2  maximum  taken  in 
isolation,  real  frequencies  would  be  limited  to  a1  <  N2QX,  and  so  k2  < 
N2ax/c2.  For  the  upper  peak,  H(tox/co='0.01rad/m,  while  for  the  lower  peak, 
the  corresponding  value  is  0.005rad/m.  The  two  families  of  curves  clearly 
display  asymptotic  approach  to  these  values.  All  curves  pass  through  the 
origin;  for  this  profile  and  source  speed,  all  modes  are  super-Froude. 

The  phase  velocities  corresponding  to  these  same  modal  curves  are 
shown  in  Figure  8.  All  phase  speeds  are  less  than  the  source  speed,  and  show 
a  monotonic  decrease  both  with  mode-  and  wave-number  in  accordance  with 
known  results7.  Figure  9  contains  plots  of  the  ratios  kXR1(ky)/ky  as  a 
function  of  ky.  For  this  combination,  it  is  evident  that  all  ratios  approach 
a  non-zero  limit  as  ky  approaches  zero.  These  real  eigenvalues  correspond  to 
the  radiating  internal  wave  system.  The  above  limiting  ratios  are  related  to 
the  angles  of  propagation  of  the  longest  horizontal  wavelengths  of  each 


Figure  7.  Plots  of  the  first  21  eigencurves  in  the  horizontal  wavenumber 
plane  for  a  source  moving  at  2.5  ms-1.  This  is  faster  than  all  of  the  phase 
speeds  of  the  internal  waves,  and  so  all  pass  through  the  origin.  Units  for 
both  axes  are  rad  nf1. 


Figure  9.  Ratios  of  kx/ky  for  the  curves  of  Figure  7,  plotted  as  a  function 
of  ky.  These  are  all  nonzero  for  ky=0,  indicating  that  all  internal  wave 
elements  have  a  cross-track  component  to  their  propagation  direction  and  so 
form  a  diverging  wake. 
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mode,  The  fact  that  they  are  all  non-zero  indicates  the  ex:  acted  results 
that  all  waves  in  the  wake  are  divergent;  fhe»-e  are  .to  transverse 
components. 


Plots  similar  to  the  above  were  also  computed  for  the  lower 
source  speed  c0  *  0.5  ms-1.  From  Figure  8,  it  is  evident  that  for  this 
source  speed,  three  sub-Froude  modes  (  c0  <  cm(G  '  would  be  expected  in 
this  case.  Plots  of  k^ky)  vs  ky  for  this  combination  are  shown  in  Figure 
10.  The  upper  three  curves  exhibit  the  expected  non-zero  kx-values  as  ky 
approaches  zero  -  these  correspond  to  modes  containing  transverse  wave 
systems  whose  longest  components  have  crests  oriented  across-track.  The 
associated  phase  velocities  for  this  source  spe^d  are  shown  in  Figure  11, 
and  show  the  effect  of  the  pole  constraint  ic0kxl  a  |kcpl  which  implies  that 
lcpl  £  1 c0 1  -  the  upper  three  curves  are  clearly  forced  to  a  limit  at  the 
source  speed  c0  *  0.5ms-1. 

(k)  Computation  of  the  Fields 

The  final  step  is  the  numerical  implementation  of  (6.14a).  As 
noted  earlier,  the  integral  itself  is  computed  by  FFT  methods,  and  so  for 
any  fixed  values  of  (x,z),  it  is  very  efficient  to  compute  many  y-values. 
However,  the  problem  of  computation  of  extensive  sets  of  field  points 
requires  some  consideration  for  efficient  computation.  Implementation  of  the 
methods  of  this  report  has  been  mainly  concerned  with  construction  of 
horizontal  (x,y)  or  vertical  (y,z)  sections  of  the  fields  at  a  fixed  value 
of  the  third  coordinate.  These  are  convenient ly  visualized  as  two- 
dimensional  "images"  on  a  video  display,  for  example. 

Since  the  output  fields  are  real,  it  suffices  to  compute  only  the 
positive  half  of  the  ky  spectra.  If  the  physical  fields  are  to  be  computed 
on  Ny  y-points  with  uniform  spacing  Ay,  the  spectra  are  then  computed 

(sampled)  r--»  ne  points  {kyJ  =  2rrj/(NyAy);  j=0 . Ny/2-l  }.  FFT  inversion 

will  then  give  a  representation  of  the  physical  fields  on  periodically 
replicated  sets  of  values  {yq,  q- 1 . . . Ny } .  For  this  set  of  ky  values,  the 


Figure  11.  Phase  speeds  for  Figure  10.  Note  that  the  moving  source 
constraint  forces  the  phase  speeds  to  be  no  greater  than  the  source  speed. 


three  sets  of  eigenvalues  (real,  imaginary,  and  zero)  corresponding  to  the 
three  classes  of  kK-poles  of  (t.15),  may  be  computed  once  for  a  given 
profile  and  source  speed,  and  saved  for  future  use. 
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This  data  set  is  independent  of  source  details  and  spatial 

coordinates.  Thus  it  may  be  used  repeatedly  for  different  source  models, 

source  depth,  and  'observation'  coordinates.  The  eigenfunctions  are  computed 

as  needed  at  run-time  -  the  amount  of  data  (  modes  x  ky-values  x  layers  )  is 

too  great  to  justify  saving.  In  general  it  will  be  desired  to  compute  the 

fields  on  a  set  of  values  {xp,  p= 1 . . . Nx]  and  {zr,  r=l...Nz}.  As  there  is 

some  numerical  expense  involved  in  construction  of  eigenfunctions,  it  is 

therefore  most  efficient  to  construct,  at  run-time,  arrays  (cf  (6.14a)  )  of 

the  general  form  {  f(x,kyj,zr);  j=0. . . Ny/2-l ,  r=l...Nz  }.  These  are  then 

inverted.  Computation  of  many  x-values  entails  computation  of  exponentials 

of  the  form  exp(-ixkx).  These  also  prove  to  be  a  significant  computational 

expense  since  there  is  one  ky-value  for  each  ku  and  mode  number  in  each  mode 

e  y 

set.  Significant  efficiency  is  therefore  realized  by  computing  fields  on  a 
uniform  set  of  x-values  with  spacing  Ax.  The  pair  of  two-dimensional  arrays 
exp[- i x1kxm(ky)  1  and  exp[ i Axkxn(ky) ]  may  be  evaluated  at  run  time,  and 
successive  x-values  obtained  by  multiplication  using 

exp[-i(x-Ax)kJ  =  exp[-ixfcj  exp[+iAxkx]. 


Sample  results  will  be  displayed  in  the  next  section. 
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8.  MAKE  EXAMPLES 


This  section  is  concerned  with  presentation  of  some  examples  of 
fluid  dynamical  fields  computed  from  the  methods  described  thus  far.  For 
purposes  of  illustration,  the  example  calculations  will  be  mostly  based  on  a 
simple  Rankine  doublet5’6,  consisting  of  a  point  source  of  strength  'm'  at 
x=+a  and  a  sink  of  equal  strength  at  x=~a.  In  the  notation  of  Equation 
(4.1),  the  mathematical  description  is 

S0(x,y)  =  m<5(y)[  5(x-a)  -  <$(x+a}],  (8.1) 

hence 

§(k)  =  (im/ir)  sin(akx). 

It  is  possible5  to  fit  the  parameters  {m,a}  to  the  diameter  and  length  of 
solid  body  with  rotationa1  symmetry  about  its  long  axis.  For  illustrative 
purposes,  the  values  of  100m  and  8m  have  been  used  for  the  length  and 
diameter  respectively.  All  examples  in  this  Section  will  be  based  on  the 
profile  of  Figure  2.  The  two  source  speeds  2.5ms-1  and  0.5ms-1  used  for 
discussion  purposes  above  will  be  considered  for  source  depths  z8  =*  35m  and 
z3  =  124m,  corresponding  to  the  centres  of  the  upper  and  lower  thermoclines. 
As  well  as  the  independent  fluid  fields  (u,v,w,£)>  the  total  energy  per  unit 
volume 


E(x.y.z)  s  Po( Z)  (  u2  +  v2  ♦  w2  +NV  )  /  2.  (8.2) 

in  the  wake  is  also  of  interest.  For  the  source  (8.1),  Equation  (6.22)  is 
appropriate.  Computations  for  all  of  the  examples  in  this  Section  are  based 
on  (6.22),  incorporating  the  first  21  terms  in  each  of  the  three  infinite 
sums  there. 


Example  1  considers  the  fields  generated  at  the  surface  by  the 
above  doublet  source  travelling  at  the  slower  speed  in  the  upper 
thermocline.  Due  to  the  boundary  conditions,  w  and  £  are  both  zero,  and  so 
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only  the  horizontal  currents  u(x,y,0)  and  v(x,y,.0)  are  present.  These  were 
computed  for  a  two-dimensional  set  of  x  and  y  values.  The  results  are 
displayed  in  Figure  12  as  pseudo-images,  wherein  the  grey  level  corresponds 
to  (signed)  current  value,  white  being  the  largest,  and  black  the  smallest. 
Computations  were  done  at  a  uniform  spacing  of  6m  in  both  x  and  y- 
directions.  There  are  512  samples  in  the  along-track  direction  and  1024 
samples  across-track.  The  physical  size  of  the  image  is  from  -2.691km  to 
0.375km  along-track,  and  from  -3.072km  to  3.066km  across-track.  The  centre 
of  the  source  is  at  x=y®0  in  the  reference  frame  chosen.  Results  have  been 
zeroed  between  the  leading  and  trailing  edges  of  the  source  in  accordance 
with  the  range  of  validity  of  Equations  (6.15).  To  produce  acceptable 
images  it  is  necessary  to  clip  the  data  to  some  desired  range,  and  map  this 
range  onto  the  available  256  grey  levels  of  the  display  unit.  In  addition, 
it  proved  desireable  to  take  the  logarithm  of  the  energy  density  image  prior 
to  the  scaling  process.  Table  I  contains  values  of  the  actual  computed 

TABLE  I. 

DATA  AND  IMAGE  RANGES  FOR  FIGURE  12. 


Field: 

u  (  mm  s-1  ) 

v  (  mm  s"1  ) 

lOQie  E  (mj  rrf3) 

minimum 

-0.408 

-1.402 

0 

maximum 

+2.690 

+1.402 

0.559 

black 

-0.120 

-0.120 

-7.0 

white 

+0.120 

+0.120 

-4.0 

maximum  and  minimum  value  for  each  of  the  fields  u,v,  and  log10E,  as  well  as 
the  values  corresponding  to  black  and  white  on  the  images,  to  which  out-of- 
range  values  have  been  cl ipped, Points  to  note  include  the  transverse  waves 
clearly  visible  in  the  along-track  velocity  component,  and  the  localized 


Figure  12.  Representations  of  the  surf<!j:e  currant  componanta  genorotad  by  c 
Rankina  doublet  source  moving  at  0.5ms  In  the  upper  therraocllne.  The  left- 
hand  Image  Is  for  the  along-track  current  and  the  central  one  for  the  cross¬ 
track  component.  Scaling  values  are  Indicated  in  Table  I.  Note  the  presence 
of  the  expected  transverse  waves  associated  with  the  first  three  modes  as 
explained  in  the  text.  The  image  size  covers  roughly  3  km  along-track  by  6 
km  cross-track  at  6m  steps  as  explained  In  the  text.  The  images  have  been 
zeroed  between  the  extremes  of  the  source,  where  the  Integration  by  residues 
breaks  down.  Note  the  localized  disturbance  near  the  source  arising  from 
Inclusion  of  the  Imaginary  modes.  The  right-hand  Image  represents  the  (  base 
10  )  logarithm  of  the  energy  density  In  the  wake. 


disturbance  In  the  vicinity  of  the  source.  Evidently,  most  of  the  energy 
density  at  the  surface  is  in  the  divergent  waves.  There  is  also  significant 
energy  density  in  the  local  disturbance. 
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Example  2  is  intended  to  illustrate  the  effect  of  source  depth  on 
the  surface  fields.  Thus  all  parameters  are  the  same  as  in  Example  1,  except 
that  the  source  is  now  travelling  in  the  lower  thermocline.  Pseudo-images  of 
the  surface  velocity  and  logie  energy  density,  analogous  to  those  of  Figure 
12,  are  presented  in  Figure  13.  Data  ranges  and  clipping  values  for  this 
image  set  are  given  in  Table  II.  The  main  point  to  note  is  the  'filtering' 
of  shorter  wave  components  by  the  interior  minimum  in  the  N2  profile.  This 
region  has  higher  evanescence  for  shorter  waves  (higher  wavenumbers),  and  so 
short  wave-length  energy  generated  in  the  lower  thermocline  is  largely 
confined  there  -  this  will  be  illustrated  further  below. 

TABLE  II. 

DATA  AND  IMAGE  RANGES  FOR  FIGURE  13. 


Field: 

u  (  mm  s-1  ) 

v  (  mm  s  1  ) 

logie  E  (mj  m-3) 

minimum 

-0.292 

-0.611 

0 

maximum 

+1.096 

+0.611 

-0.221 

black 

-0.200 

-0.060 

-7.0 

white 

+0.200 

+0.060 

-4.0 

Example  3  is  meant  to  demonstrate  the  effect  of  increasing  source 
speed,  and  uses  the  higher  value  c0  =  2.5ms-1.  Other  conditions  remain  as  in 
Example  1.  With  reference  to  Equation  (7.12),  it  is  evident  that  the  effect 
of  increasing  c0  is  to  decrease  the  parameter  sc.  Thus  the  constraint 
hyperbolae  (7.16a)  are  all  shifted  toward  smaller  values  of  p  and  hence  k  on 
account  of  (7.13).  But  since  ky  is  constant,  the  effect  is  to  reduce  kxn  at 
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FI gur«  13.  As  in  Figure  12,  but  for  a  source  moving  In  the  lower 
thermocline.  Note  the  filtering  of  shorter  waves  due  to  the  Interior  minimum 
In  the  stability  profile. 


each  ky.  Thus  every  Fourier  component  k  in  the  wake  propagates  at  an  angle 
directed  more  closely  to  the  cross  track  axis,  and  the  'V'  patterns  defined 
by  the  wake  become  narrower.  To  accommodate  this  lengthening  of  the  image, 
computations  were  done  on  a  different  mesh,  specifically  a  uniform  sample 
length  of  3m  in  both  x  and  y-direct ions .  There  are  1024  samples  in  the 
along-track  direction  and  512  samples  across-track.  The  physical  size  of  the 
image  is  from  -2.8665km  to  0.2025km  along  track,  and  from  -0.768km  to 
0.765km  across  track.  Surface  field  images  are  shown  in  Figure  14,  with  data 
range  and  clipping  values  given  in  Table  III.  There  is  clearly  no  transverse 
component  to  the  wake,  as  expected.  Also,  note  the  greatly  increased  effect 


Figura  14.  As  t n  Figure  12,  but  for  a  source  moving  at  2.5ms  .  Pixels  here 
cover  3m  by  3ra ,  while  the  total  Image  covers  roughly  3km  along-track,  and 
1.5km  ac ross-t rack .  The  localized  disturbance  Is  stronger  duo  to  the  larger 
source  strength  of  the  doublet  source.  As  expected,  only  diverging  wavos 
appear  In  the  wake. 


of  the  localized  disturbance.  The  source  strength  m  of  the  doublet  is 
roughly  proportional  to  the  source  speed,  and  therefore  results  in  an 
increased  local  flow  near  the  source. 
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TABLE  III. 

DATA  AND  I MACE  RANGES  FOR  FIGURE  14. 


Field 

u  (  mm  s_1  ) 

v  (  mm  s"1  ) 

logie  E  (mj  m~3) 

minimum 

-0.171 

-4.901 

0 

maximum 

+9.304 

+4.901 

1.636 

black 

-0.020 

-0.020 

-8.0 

white 

+0.020 

*Q.f,',fl 

-5.0 

The  next  two  examples  deal  with  vertical  sections  through  two  of 
the  above  wakes.  Both  examples  are  for  the  lower  source  speed,  and 
correspond  tc  the  two  source  depths.  The  vertical  sections  are  taken  at  x  * 
-5.5km.  The  fields  were  computed  on  the  same  y-points  as  for  Example  1,  i.e. 
1024  points  at  $m  steps.  Vertical  coordinates  consist  of  512  points  at  lm 
separation,  beginning  at  z«0.  This  vertical  range  vis  chosen  for  ease  of 
display,  and  covers  only  the  upper  portion  of  the  water  column.  For  these 
sections,  all  of  the  fluid  dynamical  fields  are  present,  and  may  be 
conveniently  presented  in  the  form  o.  images  as  above. 

Example  4  is  as  just  described,  with  the  source  in  the  lower 
tnermoclins.  The  upper  four  Images  of  Figure  15  represent  (u,v,w,£),  with 
scaling  values  given  in  Table  IV.  The  vertical  direction  represents  depth, 
while  the  horizontal  direction  represents  cross-track  distance.  These 
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TABLE  IV. 

Data  and  image  ranges  for  Figure  15. 

E1  is  for  the  deeper  source,  E2  for  the  shallow  one. 


Field 

Minimum 

Maximum 

Black 

White  . 

u  (mm/s) 

-0.107 

+0.165 

-0.107 

+0.165 

v  (mm/s) 

-0.395 

+0.395 

-0.395 

+0.395 

w  (mm/s) 

-0.641 

+0.657 

-0.641 

+0.657 

f  (cm) 

-5.977 

+6.503 

-5.977 

+6.503 

E, (mJ/s) 

1.4  x  10-6 

0.233 

1.4  x  10"6 

0.045 

E2(m«J/s) 

4.7  x  10“8 

8.545 

4.7  x  10"3 

0.489 

sections  were  computed  on  the  basis  of  1024  cross-track  points,  but  only  the 
central  512  y-points,  with  indices  from  256  to  767,  are  displayed.  Note  the 
disparate  (6-to-l)  scales  along  the  two  directions.  The  associated  energy 
density  in  the  vertical  section  is  shown  in  the  lower  left  image.  Clearly, 
the  higher  frequencies  are  confined  to  the  lower  thermocline.  There  is 
however,  significant  'leakage’  of  onger  wavelength  energy  into  the  upper 
thermocline  (which  has  a  h  gb  :r  stratification  than  the  lower  one).  This 
longer  wavelength  energy  is  ..i  fact  a'hat  makes  its  presence  felt  at  the 
surface,  as  shown  in  Figure  1.3.  Example  5,  consists  of  the  same  model,  but 
with  the  source  in  the  uoper  thermocl ire.  The  lower  right  image  of  Figure  15 
shows  the  energy  density,  ;n  the  saimj  coordinates  as  for  Example  4.  This 
represents  a  vertical  section  through  the  wake  of  Figure  12.  The  energy  is 
clearly  confined  almost  entirely  to  the  upper  tho-mocl ine.  The  velocity  and 


displacement  fields  are 
ranges  of  these  compone 


DATA  RANGES 

Field 
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it  depicted  for  Example  5,  but  for  reference,  the 
.  are  given  in  Table  V. 

TABLE  V, 

OR  THE  FIELD  COMPONENTS  OF  EXAMPLE  5 


Minimum 

Maximum 

-0.227 

0.285 

-2.76 

2.76 

-1.84 

1.74 

-15.6 

15.1 
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9.  Comparison  With  Measured  Data 

As  a  final  example,  a  comparison  with  measured  data  is 
considered.  The  details  of  the  measurements  are  discussed  in  detail 
elsewhere1.  Basically  a  surface  ship  (CFAV  Endeavour)  generated  a  divergent 
internal  wave  wake  by  steaming  in  a  straight  line  at  constant  speed.  After  a 
time,  the  ship  made  a  wide  turn,  and  made  a  measurement  cut  across  its  own 
wake.  The  data  of  interest  for  present  purposes  were  obtained  from  current 
meters  mounted  ahead  of  the  bow  at  a  nominal  depth  of  lm.  Ship  position 
information  was  obtained  from  a  trio  of  shore-based  trisponders.  The 
measured  current  velocity  data  were  reduced  to  a  component  along  the 
measurement  ship's  track.  Those  data  are  displayed  in  Figure  20  of  Reference 
1. 


Model  currents  were  computed  using  a  layered  representation  of 
the  profile  contained  in  the  above  reference,  under  all  of  the  assumptions 
and  simplifications  discussed  earlier.  The  source  type  used  was  a  more 
complicated  one  than  the  simple  Rankine  doublet.  The  measured  positions  at 
which  current  data  were  available  were  converted  into  positions  relative  to 
the  steady  state  coordinate  system  of  the  present  model.  The  synthetic 
current  data  were  then  sampled  along  this  track  and  converted  to  a  component 
along  the  measurement  track  to  produce  a  data  set  which  may  be  directly 
compared  with  the  measured  data.  The  results  of  the  comparison  are  presented 
in  Figure  16.  The  measured  data  is  the  lighter  curve,  while  the  darker  curve 
represents  the  model  data.  The  data,  in  units  of  ms-1,  are  plotted  as  a 
function  of  measurement  time.  It  is  evident  that  in  this  case  at  least, 
reasonable  agreement  was  obtained.  The  predicted  amplitudes  agree  well, 
although  the  wavelengths  are  perhaps  slightly  longer  than  measured. 
Considering  all  of  the  approximations  in  the  model,  and  the  physical 
realities  of  the  experimental  situation,  the  agreement  serves  as  at  least  an 
indication  that  linearized  Internal  wave  theory  is  capable  of  realistic 
predictions. 


1983-AUG-2 


Figure  16.  Comparison  of  internal  wave  currents  predicted  by  the  model  with 
measured  data.  The  lighter  curve  represents  data  measured  in  the  internal 
wave  wake  of  a  surface  ship.  Corresponding  data  predicted  from  the  model  are 
given  by  the  dark  curve.  The  current  is  the  component  along  the  measurement 
track  (ms-1)  as  a  function  of  measurement  time. 
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10.  Concluding  Remarks 

It  has  been  demonstrated  that  linearized  inviscid  internal  wave 
theory  is  capable  of  providing  realistic  internal  wave  field  estimates  in  an 
efficient  and  straightforward  manner.  Relatively  simple  modifications  to 
standard  normal  mode  expansion  techniques  permit  inclusion  of  both  the 
radiating  and  near-field  contributions.  These  modifications  also  provide  an 
alternative  Green's  function  representation,  which  has  a  higher  rate  of 
convergence  in  the  modal  sums.  All  pertinent  fluid  dynamical  variables  can 
be  easily  computed  together. 

It  is  also  demonstrated  that  profiles  possessing  multiple 
thermoclines  need  not  cause  major  difficulties.  For  certain  wavenumber 
regimes,  such  profiles  will  possess  regions  of  high  evanescence,  a  common 
cause  of  "instability"  in  eigenfunction  computations.  Such  profiles  will 
generally  cause  multiple  roots  to  arise  in  the  Wronskian  in  a  numerical 
scheme.  Appropriate  modifications  to  the  concept  of  an  eigenfunction 
expansion  are  readily  introduced  to  handle  such  profiles,  and  properly 
account  for  trapping  of  energy  generated  with  a  thermocline. 
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